text:
# Density of States Method at Finite Isospin Density

## 1 Introduction

The lattice QCD Monte Carlo technique has been a useful tool for clarifying the non-perturbative aspects of QCD at the zero baryon chemical potential ( $`\mu _B=0`$ ) for both finite and zero temperatures. At nonzero $`\mu _B`$, however, due to the sign problem, the standard importance sampling method fails. The Glasgow group attempted to extract information on $`\mu _B`$ at low temperature by a method based on reweighting but this method does not work at low temperatures due to the overlap problem. Recently it has been realized that at low density and finite temperature the sign problem may not be a serious numerical difficulty and a variety of approaches ( multi-reweighting, Taylor expansion, imaginary chemical potential)  have been applied to study QCD at low density and finite temperature (See for recent reviews). For larger $`\mu _B`$ we are still lacking an efficient numerical method.

In contrast to the finite $`\mu _B`$ case, QCD at isospin density has no sign problem and can be simulated by lattice Monte Carlo methods. There is an expectation that at small isospin density the system may resemble that at small baryon density. Recent studies have shown that the phase diagram at small isospin density is similar to that at small $`\mu _B`$. On the other hand, a difference between the isospin and baryon density is seen in susceptibilities of number density and derivatives of meson mass.

In lattice Monte Carlo simulations, usually importance sampling is used. This approach have been proven to be very efficient by many studies. Another approach one may take is the density of states (DOS) method. The difference between importance sampling and the DOS method is as follows. In importance sampling a simulation is performed at a fixed parameter set. To explore different parameters one must perform independent simulations111If the parameter region to be explored is narrow, then one may apply the reweighting method which does not need independent simulations.. On the other hand, in the DOS method one first determines the DOS of suitable observables and then results are obtained by performing one or a few dimensional integral. To illustrate how the DOS method works, let us consider a gauge model. The partition function for this model is given by $`Z={\displaystyle [dU]\mathrm{exp}(\beta S_g[U])},`$ where $`S_g[U]`$ is a gauge action and $`\beta `$ is a coupling constant. If we define the DOS by $`n(E)={\displaystyle [dU]\delta (ES_g[U])},`$ then the partition function is rewritten as $`Z={\displaystyle 𝑑En(E)\mathrm{exp}(\beta E)}.`$ If we consider the average value of the gauge action, it is given by $`<E>=𝑑En(E)E\mathrm{exp}(\beta E)/Z`$. This is a one-dimensional integral of $`E`$ and once we obtain $`n(E)`$, we can evaluate this easily at $`any`$ $`\beta `$. No extra simulation is needed to calculate it at various $`\beta `$, which is considered to be an advantage of this method.

The DOS method has been applied for gauge theories: Z(n), SU(2) and SU(3) models. Although the inclusion of dynamical fermions is computationally difficult, in Ref. the DOS method was applied for QED, where the method is called the microcanonical fermionic average method. Recently, Luo extended the idea of the DOS method to QCD and emphasized that once the eigenvalues of the Dirac operator and the DOS of the plaquette energy are determined, one can evaluate the thermodynamic quantities derived from the partition function at $`any`$ quark mass and flavor.

In this study we apply the DOS method for QCD at isospin density on a $`4^4`$ lattice and demonstrate how the DOS method works. Since the isospin system has no sign problem, results from the DOS method can be compared with those from the standard method. In Sec.2 we give general formulas for forming the DOS including dynamical fermions. In Sec.3 we give simulation details. Results are presented in Sec.4. We summarize our results in Sec.5.

## 2 Density of States Method

### 2.1 General Formulas

In lattice QCD Monte Carlo simulations, usually we aim at obtaining the average values of observables:

$$O=\frac{1}{Z}[dU]O[U]det\mathrm{\Delta }(m_q,\mu )^{N_f/4}\mathrm{exp}(\beta S_g[U]),$$
(1)
where $`N_f`$ is the number of flavors and $`\mathrm{\Delta }(m_q,\mu )`$ is assumed to be a staggered fermion matrix at quark mass $`m_q`$ and at chemical potential $`\mu `$:

$`\mathrm{\Delta }(m_q,\mu )_{i,j}=m_q\delta _{i,j}`$ $`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\nu =\widehat{1},\widehat{2},\widehat{3}}{}}\eta _{i,\nu }(U_{i,\nu }\delta _{i,j\nu }U_{i\nu ,\nu }^{}\delta _{i,j+\nu })`$
$`+`$ $`{\displaystyle \frac{1}{2}}\eta _{i,\widehat{4}}(e^\mu U_{i,\widehat{4}}\delta _{i,j\widehat{4}}e^\mu U_{i\widehat{4},\widehat{4}}^{}\delta _{i,j+\widehat{4}}).`$
$`S_g[U]`$ is the standard Wilson gauge action,

$$S_g[U]=\frac{1}{N_c}\underset{p}{}ReTr(U_p),$$
(3)
where $`U_p`$ is the plaquette and $`N_c=3`$ for $`SU(3)`$ gauge theory. For Wilson fermions, $`N_f/4`$ in eq.(1) should be replaced with $`N_f`$. The standard means of dealing with this multi-dimensional integral is to do the importance sampling, i.e. configurations are generated with a measure $`dUdet\mathrm{\Delta }^{N_f/4}e^{\beta S_g}`$ and $`O`$ is given by an average over the configurations as $`O\frac{1}{N}_{i=1}^NO[U_i]`$, where $`N`$ is the number of configurations.

The DOS method rewrites the partition function and reduces it to a few dimensional integral. Let us define the DOS with parameters $`E_i`$ ($`i=1,\mathrm{},k`$) as

$$n(E_1,E_2,\mathrm{},E_k)=[dU]g(U)\mathrm{\Pi }_i\delta (E_ix_i(U)),$$
(4)
where $`x_i(U)`$ is an operator associated with $`E_i`$ and $`g(U)`$ is introduced to generalize the DOS further with $`g(U)1`$. Using the DOS, we obtain

$$O=\frac{1}{Z_n}[\mathrm{\Pi }_idE_i]n(E_1,E_2,\mathrm{},E_k)Odet\mathrm{\Delta }^{N_f/4}e^{\beta S_g}/g(U)_{E_1,E_1,\mathrm{},E_k},$$
(5)
where

$$Z_n=[\mathrm{\Pi }_idE_i]n(E_1,E_2,\mathrm{},E_k)det\mathrm{\Delta }^{N_f/4}e^{\beta S_g}/g(U)_{E_1,E_2,\mathrm{},E_k}.$$
(6)
Here, $`A_{E_1,E_2,\mathrm{},E_k}`$ stands for an average of $`A`$ on configurations generated with the measure $`[dU]g(U)\mathrm{\Pi }_i\delta (E_ix_i(U))`$, or one can write it as

$$A_{E_1,E_2,\mathrm{},E_k}=\frac{1}{n(E_1,E_2,\mathrm{},E_k)}𝑑Ug(U)\mathrm{\Pi }_i\delta (E_ix_i(U))A(U).$$
(7)
For $`g(U)=1`$, $`A_{E_1,E_2,\mathrm{},E_k}`$ becomes a microcanonical average at fixed $`E_1,E_2,\mathrm{},E_k`$.

### 2.2 $`g(U)=1`$

Here we give the formulas used in the present study. We consider the case of one parameter and choose $`x(U)=S_g[U]`$. Although there are many possibilities to choose $`x(U)`$, this choice with $`g(U)=1`$ is the one used in Ref. and is useful for our purpose. Setting $`g(U)=1`$, we obtain

$$n(E)=[dU]\delta (6VES_g[U]),$$
(8)
$$O=\frac{1}{Z_n}𝑑En(E)e^{\beta 6VE}Odet\mathrm{\Delta }(\mu )^{N_f/4}_E,$$
(9)
and

$$Z_n=𝑑En(E)e^{\beta 6VE}det\mathrm{\Delta }(\mu )^{N_f/4}_E,$$
(10)
where $`V`$ is the number of lattice sites and $`E`$ is the plaquette energy. These are the basic formulas used in our study. Our task includes three numerical calculations: (i)$`n(E)`$, (ii) $`det\mathrm{\Delta }(\mu )^{N_f/4}_E`$ and (iii) $`Odet\mathrm{\Delta }(\mu )^{N_f/4}_E`$. Luo  argued that if one stores the eigenvalues of the fermion matrix for all configurations, then one can evaluate $`det\mathrm{\Delta }(\mu )^{N_f/4}_E`$ for any quark mass and flavor. Let $`\lambda _i(\mu )`$ be the i-th eigenvalue of the massless fermion matrix $`\mathrm{\Delta }(m_q=0)`$ on a configuration with $`E`$. Then we obtain

$$det\mathrm{\Delta }(\mu )^{N_f/4}_E=(\underset{i}{\overset{N_cV}{}}(\lambda _i(\mu )+m_q))^{N_f/4}_E.$$
(11)
Since (iii) contains $`O`$, in general it is not calculable for any quark mass and flavor. However the chiral condensate $`\overline{\psi }\psi `$ which is obtained as the trace of the inverse fermion matrix can also be given with the eigenvalues. Namely, for $`\overline{\psi }\psi `$ we obtain

$$(\frac{1}{V}Tr\mathrm{\Delta }(\mu )^1)det\mathrm{\Delta }(\mu )^{N_f/4}_E=\frac{1}{V}\underset{i}{}\frac{1}{\lambda _i(\mu )+m_q}(\underset{i}{\overset{N_cV}{}}(\lambda _i(\mu )+m_q))^{N_f/4}_E.$$
(12)
### 2.3 $`g(U)1`$

Depending on the purpose, one may take $`g(U)1`$. In Ref. the DOS for the complex phase of the fermion determinant was calculated. The DOS was defined as

$$n(E)=[dU]g(U)\delta (E\theta (U)),$$
(13)
where $`g(U)`$ is set to $`|det\mathrm{\Delta }|e^{\beta S_g(U)}`$, and $`\theta (U)`$ is the complex phase of the determinant. Then the expectation value of $`O`$ is given by

$$O=_\pi ^\pi 𝑑EO_En(E)e^{iE}/_\pi ^\pi 𝑑En(E)e^{iE}.$$
(14)
In Ref. the DOS for the number density was used and $`g(U)`$ is also set to $`|det\mathrm{\Delta }|e^{\beta S_g(U)}`$.

In these definitions parameters $`\beta `$ and $`m_q`$ are absorbed in $`n(E)`$ and we cannot vary $`\beta `$ and $`m_q`$. Therefore, in this study we do not use these definitions with $`g(U)1`$.

## 3 Simulations

For the case of $`N_f=2`$, we have two chemical potentials $`\mu _u`$ and $`\mu _d`$ associated with the $`u`$ and $`d`$ quarks respectively. If we introduce the isospin chemical potential $`\mu _I`$ as $`\mu _I\mu _u=\mu _d`$, we obtain

$$[det\mathrm{\Delta }(\mu _u)\mathrm{\Delta }(\mu _d)]^{1/4}=[det\mathrm{\Delta }(\mu _I)\mathrm{\Delta }(\mu _I)]^{1/4}=|det\mathrm{\Delta }(\mu _I)|^{1/2},$$
(15)
which is positive in the Monte Carlo measure. Therefore, for this $`\mu _I`$ the sign problem is absent and one can perform Monte Carlo simulations. The natural generalization of the isospin chemical potential to $`N_f2`$ is that one takes $`det\mathrm{\Delta }(\mu _I)^{N_f/4}|det\mathrm{\Delta }(\mu _I)|^{N_f/4}`$.

We follow the implementation developed in Ref.. The DOS $`n(E)`$ in eq.(8) can be obtained using the quenched data as

$$\frac{\mathrm{ln}n(E)}{V}=6_0^E𝑑E^{}\beta (E^{})+const.$$
(16)
First we make quenched simulations on a $`4^4`$ lattice and determine the coupling constant $`\beta (E)`$ as a function of the plaquette $`E`$ ( Figure 1 ). Then we integrate $`\beta (E)`$, numerically interpolating the data according to the trapezoidal rule. Figure 2 shows the result of $`\mathrm{ln}n(E)/V`$.

The most time consuming part of our method is the calculations of the microcanonical averages $`det\mathrm{\Delta }(\mu _I)^{N_f/4}_E`$ and $`Odet\mathrm{\Delta }(\mu _I)^{N_f/4}_E`$ as a function of $`E`$, which contain the eigenvalue calculations. In order to generate configurations at $`E`$ we use the over-relaxation method. Starting from a configuration with $`E`$, we have made the over-relaxation update and saved 100 configurations at each $`E`$. Each configuration is separated by 100 over-relaxation updates. About 30 values of $`E`$ are chosen in $`E[0.05,0.95]`$. For each configuration we calculate the eigenvalues of the fermion matrix and store them. Using these eigenvalues we can evaluate $`det\mathrm{\Delta }(\mu _I)^{N_f/4}_E`$ as in eq.(11) for any quark mass and flavor. Examples of $`det\mathrm{\Delta }(\mu _I)^{N_f/4}_E`$ are shown in Figure 3. Eq.(12) is evaluated similarly. These data are interpolated by polynomials when we perform eqs.(8)-(10).

## 4 Results

Figures 4 and 5 compare results of $`\overline{\psi }\psi `$ between the DOS method and the R-algorithm. We see a good agreement between them in a wide range of $`\beta `$. A small difference is seen in the phase transition region where the $`\overline{\psi }\psi `$ changes rapidly. In such a region we probably need careful analysis.

Figure 6 shows $`\overline{\psi }\psi `$ for different $`N_f`$ at $`m_q=0.025`$ and $`\mu _I=0.2`$(left), and at $`m_q=0.05`$ and $`\mu _I=0.25`$(right) as functions of $`\beta `$. One can see that the critical coupling where the phase transition occurs decreases as $`N_f`$ increases.

Figure 7 shows $`\overline{\psi }\psi `$ for various $`\mu _I`$. In the low-temperature phase $`\overline{\psi }\psi `$ decreases as $`\mu _I`$ increases. On the other hand, in the high-temperature phase no visible difference can be seen among $`\overline{\psi }\psi `$. Since the DOS method can explore various parameter space easily it is considered to be useful for exploring a wide parameter space and for seeing a rough phase diagram.

In the DOS method we can take various combinations of parameters. Let us consider the case of $`N_f=1+1`$ with non-degenerate quark masses $`m_1`$ and $`m_2`$. In this case we must calculate the following microcanonical averages:

$$|det\mathrm{\Delta }(m_1)|^{N_f/4}|det\mathrm{\Delta }(m_2)|^{N_f/4}_E,$$
(17)
$$\overline{\psi }\psi (m_{i=1,2})|det\mathrm{\Delta }(m_1)|^{N_f/4}|det\mathrm{\Delta }(m_2)|^{N_f/4}_E.$$
(18)
Since the eigenvalues are stored, it is easy to calculate these microcanonical averages. On the other hand, in the conventional algorithm such as R-algorithm one needs a differently implemented program to simulate $`N_f=1+1`$. Such a program can be implemented but may become intricate.

Figure 8 shows $`\overline{\psi }\psi `$ for $`N_f=1+1`$ with different quark masses ( $`m_1=0.05`$ and $`m_2=0.025`$ ) and at $`\mu _I=0.2`$. Since the two quarks have different masses, we plot two $`\overline{\psi }\psi (m_q)`$ for $`m_1`$ and $`m_2`$. Similarly, we can also consider non-degenerate isospin chemical potentials ( $`\mu _1`$ and $`\mu _2`$ ). In this case we calculate

$$|det\mathrm{\Delta }(\mu _1)|^{N_f/4}|det\mathrm{\Delta }(\mu _2)|^{N_f/4}_E,$$
(19)
$$\overline{\psi }\psi (\mu _{i=1,2})|det\mathrm{\Delta }(\mu _1)|^{N_f/4}|det\mathrm{\Delta }(\mu _2)|^{N_f/4}_E.$$
(20)
Figure 9 shows $`\overline{\psi }\psi `$ for $`N_f=1+1`$ with different chemical potentials ( $`\mu _1=0.2`$ and $`\mu _2=0.3`$ ) at $`m_q=0.025`$.

Next we show the results of the number density $`n_d`$ defined by

$$n_d=\frac{1}{V}\frac{\mathrm{ln}Z}{\mu }=\frac{1}{V}Tr\frac{1}{\mathrm{\Delta }(\mu )}\frac{\mathrm{\Delta }(\mu )}{\mu }.$$
(21)
Since this can not be evaluated with the eigenvalues of $`\mathrm{\Delta }(\mu )`$ only, in order to obtain $`O|det\mathrm{\Delta }(\mu )|^{N_f/4}_E`$ for the number density we must calculate the number density on each configuration.

Usually Eq.(21) is evaluated by the noise method. Using noise vectors $`R_i`$ having the property $`R_i^{}R_j=\delta _{ij}`$, the trace calculation can be replaced by

$$Tr\frac{1}{\mathrm{\Delta }(\mu )}\frac{\mathrm{\Delta }(\mu )}{\mu }=\frac{1}{N}\underset{i=1}{\overset{N}{}}R_i^{}\frac{1}{\mathrm{\Delta }(\mu )}\frac{\mathrm{\Delta }(\mu )}{\mu }R_i,$$
(22)
where $`N`$ is the number of noise vectors. Although this is true for $`N\mathrm{}`$, we find that the convergence of the noise method is extremely slow on a $`4^4`$ lattice. A similar result was reported in Ref.. Figure 10 shows the convergence measured as relative error as a function of the number of $`Z_2`$ noise vectors. Typically, $`O(1000)`$ noise vectors are needed to have a reasonable value. In this analysis, instead of using the noise method we calculated the number density exactly by calculating each diagonal element of $`{\displaystyle \frac{1}{\mathrm{\Delta }(\mu )}}{\displaystyle \frac{\mathrm{\Delta }(\mu )}{\mu }}`$. In general such calculations are computationally costly and should be avoided. However, our lattice size is sufficiently small to perform the exact calculation. Thus, here we have adopted the exact calculation.

Figure 11 shows the number density as a function of $`\mu _I`$. The results from the R-algorithm are also plotted in the figure. The two results are in good agreement.

## 5 Summary and outlook

We have given the general formulas of the DOS method including dynamical fermions. The case of ”$`g(U)=1`$” corresponds to that used by Luo. Based on the implementation by Luo, we have calculated $`\overline{\psi }\psi `$ and the number density at finite isospin densities by the DOS method and made a comparison with results from the R-algorithm. The two results were found to be in good agreement. The DOS method can be applied for various combinations of parameters. We have calculated $`\overline{\psi }\psi `$ for various values of $`N_f`$ and $`\mu _I`$, and also for non-degenerate quark masses and for different isospin chemical potentials. Especially it is emphasized that for non-degenerate quark masses and for different isospin chemical potentials it is not easy to perform Monte Carlo simulations by the conventional algorithm as the R-algorithm but these calculations are easily performed in the DOS method by keeping all eigenvalues of the fermion matrix.

The limitation of the DOS method may appear on a large lattice. The measurement of the DOS method is done in the microcanonical ensemble and $`det\mathrm{\Delta }`$ is treated as an observable. Since the microcanonical and full ensembles are very different, $`det\mathrm{\Delta }`$ is expected to fluctuate largely as the volume of the system increases, which may limit the available lattice to a small one.

In principle, the DOS method can be applied for QCD with a baryon chemical potential. However, there still remains the sign problem. For example, for $`N_f=4`$ we have $`det\mathrm{\Delta }(\mu _B)_E=|det\mathrm{\Delta }(\mu _B)|e^{i\theta }_E`$ and if $`e^{i\theta }`$ fluctuates significantly, which is expected to occur for larger $`\mu _B`$, one cannot obtain meaningful values for the microcanonical average. We attempted to calculate $`|det\mathrm{\Delta }(\mu _B)|e^{i\theta }_E`$ but for $`\mu _B>0.2`$ we could not obtain statistically meaningful values222Here $`\mu _B`$ is defined as the quark chemical potential.. On the other hand, for $`\mu _B<0.2`$ the results were stable but there is no visible difference between the results of $`\overline{\psi }\psi `$ at $`\mu _I`$ and at $`\mu _B`$; for $`\mu _B<0.2`$ the effect of the complex phase is small and thus we do not see any difference. This is consistent with the results of Refs..

In this study for each value of $`\mu `$ we calculated the eigenvalues. However, if we use a matrix reduced to two time slices, we can calculate the determinant at any $`\mu `$. The determinant for this case is expressed as

$$det\mathrm{\Delta }(\mu )=det(P+e^{N_t\mu })\times e^{3V\mu },$$
(23)
where $`P`$ is a $`\mu `$-independent matrix and $`N_t`$ is the number of lattice sites in the time direction. Keeping the eigenvalues of $`P`$, we can calculate the determinant at any $`\mu `$. In this case, however, the tradeoff is that the determinant is not calculable at any $`m_q`$. Thus, if one wishes to obtain the behavior with $`\mu `$ at fixed $`m_q`$ this reduction method may be useful.

## Acknowledgments

These calculations were carried out on the SX-5 at RCNP, Osaka University. The author would like to thank Ph. de Forcrand, X.Q.Luo and A.Nakamura for helpful comments and discussions.


format: text

text:
# Chemical compositions of four barium stars

## 1 Introduction

The “barium stars”, or “Ba ii stars” as they were first called, were defined as a class of peculiar red giants by Bidelman & Keenan (1951). These G and K giants exhibit enhanced features of Ba ii, Sr ii, CH, CN, and sometimes C2 lines. Many qualitative studies on Ba ii stars have been developed (e.g. Garstang 1952; Burbidge & Burbidge 1957; Danziger 1965; Pilachowski 1977; Tomkin & Lambert 1983; Smith 1984; Kovacs 1985; Za$`\stackrel{ˇ}{c}`$s 1994). The results confirmed that, in addition to Ba and Sr, other heavy elements are also enhanced, e.g. Y, Zr, La, Ce, Pr, Nd and Sm.

Some researchers suggested that the overabundances of Ba stars were produced in the interior of the Ba ii star, and subsequently mixed to the surface (e.g. Burbidge & Burbidge 1957 for HD 46407). However, the absence of the unstable nucleus 99Tc ($`\tau _{\frac{1}{2}}=2\times 10^5`$ yrs) supply evidence that they are not presently undergoing nucleosynthesis activity and the third dredge-up process like asymptotic giant branch (AGB) stars.

After Burbidge et al. (1957) suggested that elements heavier than iron are synthesized in the interior of AGB stars, it is generally believed that Ba stars belong to binary systems and their heavy-element overabundances are produced by accreting the matter ejected by the companions (the former AGB stars, now evolved into white dwarfs). Many researchers have studied the binarity or heavy-element abundances of Ba stars (McClure et al. 1980; McClure 1983; McClure & Woodsworth 1990; Boffin & Jorissen 1988; Jorissen & Mayor 1992; Jorissen et al. 1998; Liang et al. 2000). The mass exchange took place about $`1\times 10^6`$ years ago, so the 99Tc produced in the original thermal pulse (TP) AGB stars have decayed. The accretion may be in the forms of wind accretion, disk accretion or common envelope ejection (Han et al. 1995; Jorissen et al. 1998; Liang et al. 2000).

Therefore, it is very important to combine abundance analysis and binary characters of Ba stars to study their properties and formation scenarios.

At present, there is a large sample of Ba stars with measurements of orbital elements (Carquillat et al. 1998; Udry et al. 1998a, 1998b; Jorissen et al. 1998), absolute magnitudes and kinematics (Gómez et al. 1997; Mennessier et al. 1997). However, the corresponding heavy-element abundances have not been obtained from observations. To explore their formation scenario, it is necessary to study their observed abundances, combined with their orbital elements.

Moreover, using our angular momentum conservation model of wind accretion of Ba stars formation, and combining the AGB stars nucleosynthesis (Liang et al. 2000; Liu et al. 2000), we can calculate the theoretical heavy-element abundances of Ba stars. Thus we can understand the observed abundance patterns of the sample stars from theory.

In addition, with the high precision Hipparcos data, the precise photometric parameters, improved methods to determine stellar atmospheric parameters and developed stellar evolutionary tracks study, it is possible to allow to obtain the masses and atmospheric parameters of these stars. For Ba stars, combining masses, element abundances and orbital parameters provides a good understanding of their properties and formation. In this paper, we discuss the properties of four Ba stars by combining these parameters.

This paper is organized as follows. In Sect. 2–4, we describe the observations and analysis methods in details, and present the derived atmosphere parameters and masses. Detailed spectral lines analysis and equivalent widths (EWs) are also presented. In Sect. 5, we show the element abundance patterns and analyze the abundance uncertainties. In Sect. 6, we calculate the theoretical abundance trends of s-process elements in Ba stars, and compare them with the observed patterns of our sample stars. Detailed discussions are given in Sect. 7. In Sect. 8, we summarize this paper and describe possible future studies.

## 2 Observations and data reduction

The high resolution spectra of four barium stars were obtained with the Coudé Echelle Spectrograph and a 1024$`\times `$1024 TeK CCD attached to the 2.16m telescope at the National Astronomical Observatories (Xinglong, China) in 2001. The red arm of the spectrograph with a 31.6 grooves/mm grating was used. With a 0.5 mm slit, the resolving power is of the order of 40 000 in the middle focus camera system. The total wavelength coverage was 5800$``$9000 Å . The signal-to-noise (S/N) ratio is greater than 150 over the whole region. A detailed description of technical aspects of the spectrograph can be found in Zhao & Li (2001).

The data reductions were carried out with a standard MIDAS package. The procedure is: order identification, background subtraction, flat-field correction, order extraction, wavelength calibration, radial velocity shift correction, and continuum normalization. Bias, dark current and scattered–light corrections are included in the background subtraction. The pixel-to-pixel sensitivity variations were corrected by using the flat-field. The wavelength calibration was based on thorium-argon lamp spectra. In Fig. 1, we present a portion of spectrum of the sample star HD 98839 in the region near Ba ii $`\lambda `$6141.727 line.

The measurements of spectral line EWs were done by applying two different methods: direct integration of the line profile and Gaussian fitting. The latter is preferable in the case of faint lines, but unsuitable for the strong lines in which the damping wings contribute significantly to the equivalent width. The final EWs are weighted averages of these two measurements, depending on the line intensity (see Zhao et al. 2000 for details).

## 3 Stellar atmospheric parameters

Model atmospheres are specified by four parameters: effective temperature, surface gravity, overall metallicity and microturbulent velocity.

Effective temperature $`T_{\mathrm{eff}}`$ is determined from the $`uvby`$ system index ($`by`$) and \[Fe/H\] using the empirical calibration of Alonso et al. (1999). This calibration is based on a large sample of field and globular cluster stars roughly covering spectral type from F0 to K5, and derived by applying the infrared flux method. The photometric data of sample stars are taken from Olsen (1993) and Hauck & Mermilliod (1998). The error of the photometric data is $`\sigma (by)`$=0.004 mag. Adopting the mean $`\sigma `$(Fe/H)=0.14 dex from the spectroscopic analysis, the statistical error of $`T_{\mathrm{eff}}`$ is estimated to about $`\pm `$50 K. Considering a possible error of $`\pm `$80 K in the calibration, and other factors, we estimate the error in $`T_{\mathrm{eff}}`$ as high as 100 K for our sample stars.

The previous method to determine gravity log $`g`$ requires that Fe i and Fe ii lines give the same iron abundance. But it is well known that the derivation of iron abundance from Fe i and Fe ii lines may be affected by many factors such as unreliable oscillator strengths, possible non $`localthermodynamicequilibrium`$ (non-LTE) effects and uncertainties in the temperature structure of the model atmospheres. From Hipparcos parallaxes, we can determine more reliable gravities. Using the relations:

$$log\frac{g}{g_{}}=log\frac{M}{M_{}}+4log\frac{T_{\mathrm{eff}}}{T_{\mathrm{eff},}}+0.4(M_{\mathrm{bol}}M_{\mathrm{bol},})$$

and

$$M_{\mathrm{bol}}=V+5+5log\pi +BC,$$

where $`M`$ is the stellar mass, $`M_{\mathrm{bol}}`$ the absolute bolometric magnitude, $`V`$ the visual magnitude, $`BC`$ the bolometric correction, and $`\pi `$ the parallax. We adopt solar value log $`g_{}`$=4.44, $`T_{\mathrm{eff},}`$=5770 K, $`M_{\mathrm{bol},}`$=4.75 mag. The parallax $`\pi `$ is taken from the Hipparcos Satellite observations (ESA 1997). Stellar mass was determined from the position of the star in $`M_{\mathrm{bol}}`$log $`T_{\mathrm{eff}}`$ diagram. The corresponding errors were estimated by changing $`T_{\mathrm{eff}}`$, \[Fe/H\] and $`M_{\mathrm{bol}}`$ according to its derivations (more details can be found in Chen et al. 2000). We adopt the stellar evolution tracks given by Girardi et al. (2000). The bolometric correction, $`BC`$, is determined using the empirical calibration of Alonso et al. (1999). The uncertainty of log $`g`$ determined by this method mainly comes from the uncertainties in the stellar distance, stellar mass, effective temperature and bolometric correction. We estimate the error in log $`g`$ to be in the range of 0.1$``$0.3 dex for our sample stars (Table 2 and the Col. 4 in Table 5).

The initial metallicities of the program stars were taken from the literature if available. Otherwise, we obtain an initial value by judging from the spectra and the color index. The final model metallicity is derived from consistency with the other parameters in the abundance calculation.

Microturbulent velocity $`\xi _t`$ was determined from the abundance analysis by requiring a zero slope of \[Fe/H\] vs. EWs. The values of the sample stars are given in Table 2. The error of $`\xi _t`$ is about 0.3 km s-1.

Some characteristic parameters of four barium stars are presented in Table 1. The columns are: star name, spectral type, $`V`$ magnitude, ($`BV`$), ($`by`$), parallax $`\pi `$, uncertainty of parallax, bolometric correction and stellar mass. The atmospheric parameters of the sample stars are summarized in Table 2 including star name, effective temperature, surface gravity, microturbulent velocity and metallicity.

## 4 Model atmospheres and spectral lines

We adopt the model atmospheres generated by ATLAS9 code (Kurucz 1993) to do the abundance analysis: these are LTE, plane-parallel, line-blanketed models. Abundances are determined by using the input atmospheric parameters given in Table 2 and the measured EWs. We usually choose the lines with EWs$`=`$20$``$120 mÅ, because the weaker lines would increase random error and possibly some systematic overestimates, while the stronger lines are very sensitive to the microturbulence and damping. We use the updated oscillator strengths of spectrum lines taken from the NIST database (http://physics.nist.gov). But the log $`gf`$ values of n-capture process elements are taken from the related references (Hannaford et al. 1982; Bi$`\stackrel{´}{e}`$mont et al. 1981; Weise & Martin 1980; Luck & Bond 1991; Bi$`\stackrel{´}{e}`$mont et al. 1982). All of the lines used in determining abundances are given in Table 3, including line and wavelength, excitation potential $`\chi `$, oscillator strengths log $`gf`$ and EWs. For Ba ii $`\lambda `$5853.688 line in HD 98839, we find its EW is 171 mÅ, which is much lower than 259 mÅ obtained by Fernández-Villacañas et al. (1990). The EW values of Ba ii $`\lambda \lambda `$6141.727, 6496.908 lines in this star are higher than 200 mÅ obtained in this study.

## 5 Abundances and their uncertainties

### 5.1 Abundance results

The final abundances log $`ϵ`$ (in the usual scale log$`X`$(H) = 12.0) and the corresponding induced \[$`X`$/Fe\] values are presented in Table 4, where N refers to the line number used in calculation. The solar abundances are taken from Grevesse & Sauval (1998).

Fig. 2(a)-(d) shows the detailed abundance results of the sample stars. The corresponding errors are taken from Col. 6, $`\sigma _{tot}`$, in Table 5. Detailed analysis about errors is discussed in next section (Sect. 5.2).

### 5.2 Errors in resulting abundances

Uncertainties of the abundances mainly come from equivalent widths measurement, log $`gf`$ values and stellar atmospheric parameters.

Uncertainties in the EWs are set essentially by the S/N and the resolution of the spectra. An expression for the intrinsic accuracy of an equivalent width is provided by Cayrel (1988). Applied to our spectra, having R=40,000 and typical S/N of 150, we expect uncertainties in the EWs of about 2 mÅ. We also consider the errors obtained by measuring EWs of spectral lines several times in the total uncertainties. For an element represented by N lines, the error is decreased by a factor $`\sqrt{N}`$. Thus the errors in element abundances from EWs are calculated and illustrated on Col. 2 in Table 5, $`\sigma _{EW}/\sqrt{N}`$.

The errors caused by atmospheric parameters can be obtained through changing the parameters $`T_{\mathrm{eff}}`$, log $`g`$, $`\xi _t`$ individually. Column 3, 4, 5 in Table 5 show the effects on the derived abundances by changing +100 K in $`T_{\mathrm{eff}}`$, +0.3, +0.2 or +0.1 dex in log $`g`$, +0.3 km s-1 in $`\xi _t`$.

Model metallicity uncertainties scarcely affect the element abundances, and they can be ignored.

The uncertainty in $`gf`$-values of the Fe i and Fe ii lines is about 0.05 dex. For the elements up to Ni, the uncertainties are about 0.02$``$0.10 dex, and for the heavier elements, the errors are $``$0.20 dex (Smith et al. 1995; Lambert et al. 1996; Junqueira & Pereira 2001).

The combined abundance errors caused by the above-mentioned uncertainties in $`gf`$-values, EWs-values and stellar atmospheric parameters, are the total errors in abundances, which are given in the last columns in Table 5.(1)-(4) by $`\sigma _{tot}`$.

### 5.3 Abundance pattern analysis

Among the four sample stars, HD 98839 shows the highest metallicity with \[Fe/H\]=+0.16, which is slightly higher than solar metallicity. HD 50082 shows the lowest metallicity with \[Fe/H\]=–0.40. The other two stars (\[Fe/H\] is –0.06 for HD 27271 and –0.28 for HD26886) are nearer to solar metallicity than HD 98839 and HD 50082. However, it is clear that all of them are disk stars.

The \[$`\alpha `$/Fe\] of sample stars are similar to solar. And it is clear that \[$`\alpha `$/Fe\] patterns are lower than those of the two halo Ba stars (\[Fe/H\]=–1.72, –1.43) studied by Junqueira & Pereira (2001). Their two halo Ba stars show the overabundant \[$`\alpha `$/Fe\] relative to solar. The reason is that massive stars produce these $`\alpha `$ elements, which have been ejected into interstellar medium with the Type II supernova explosion of the massive stars since the early stage of our Galaxy. Most of element Fe is produced through Type Ia supernova explosion, which occurs in the close binary system composed by the longer lifetime, intermediate- and low-mass stars (Timmes et al. 1995; Samland 1998; Liang et al. 2001). Thus, there is abundant element Fe in our Galaxy at its present stage. The \[$`\alpha `$/Fe\] patterns of our sample Ba stars also show that all of these four stars are disk stars. In addition, it should be noticed that the \[O/Fe\] values of our three samples will be increased by 0.14 when we use the recent solar O abundance log $`ϵ`$(O)=8.69 obtained by Allende Prieto et al. (2001) based on \[O i\] $`\lambda `$6300 forbidden line. That is, the \[O/Fe\] ratio in HD 26886, HD 27271 and HD 98839 are 0.50, 0.53 and 0.29, respectively. But this increase value is still smaller than the corresponding abundance errors in the three sample stars: 0.22, 0.18 and 0.21, respectively.

For Na and Al, these two odd elements, their abundances are similar to Mg abundance, because these nuclei are thought to be synthesized in the same process, namely, carbon burning in massive stars.

The iron group elements follow Fe closely though with some scatter. Generally, Cr and Mn are slightly overabundant, and Ti, V are slightly underabundant. Especially, Ni, the closest element to Fe, is very similar to the Fe abundance in the sample stars.

However, it is not clear whether the larger scatter of $`\alpha `$ elements and iron elements in HD 50082 is the natural results of its strong Ba property or not.

The most important results are the obvious overabundance of s-process elements ($`Z>`$56), Y, Zr and Ba, La, Eu (two peaks), in all of the four Ba stars (without appropriate Y and Zr lines are used in HD 98839 due to spectral qualities). Especially, the second peak s-process elements, Ba and La, are in the highest abundances. The origin of the two peaks is the result of neutron magic number 50 and 82 nuclei on the path of s-process nucleosynthesis occurred in the interiors of AGB stars (Clayton et al. 1961; Seeger et al. 1965), the progenitor companions of these Ba stars. There may be a third peak on element Pb corresponding to the neutron magic number 126. But it needs high quality blue spectra to get the Pb i $`\lambda `$$`\lambda `$3683, 4057 lines to obtain the Pb abundance.

To further degree, in the three sample Ba stars, the first peak s-elements, Y and Zr, show lower abundances than the second peak s-elements Ba, La etc. This trend also was shown in some other Ba stars given by Za$`\stackrel{ˇ}{c}`$s (1994). This \[hs/ls\]$`>`$0 trend (“hs” refers to the second peak elements, “ls” refers to the first peak elements) shows that the neutron exposures experienced in their AGB progenitor companions are strong enough to produce the more abundant “hs” elements than the “ls” elements. If we plot our Ba samples and the corresponding AGB progenitor companions on Fig. 1 of Liang et al. 2000), or Fig. 6 of Busso et al. (1995) (\[ls/Fe\] vs. \[hs/ls\]), these stars should be in the \[hs/ls\]$`>`$0 regions in the two figures. Our theoretical calculations given in next section certify this point.

This discrepancy between the two peak elements may be caused by the different neutron exposures in the companion AGB stars. Namely, the higher neutron exposure in the interior of AGB stars benefits to the production of the heavier s-process elements (Liang et al. 2000; Busso et al. 1995).

Comparing the abundance patterns of the two “typical” Ba star HD 27271 and HD 50082, we notice that the differences between the two peak elements (Y, Zr and Ba, La) are higher in HD 50082 than in HD 27271. It may be that the companion AGB star of HD 50082 experiences stronger neutron exposure in s-process nucleosynthesis. As a result, the abundances of all of the s-process elements of HD 50082 are higher than those of HD 27271. Our theoretical calculations certify this analysis (see Sect. 6).

As for calling HD 50082 and HD 27271 as “typical” Ba stars, the main reason is their masses. According to the scale defined by Warner (1965), there are three mild Ba stars in our samples: HD 26886, HD 27271 and HD 98839, and one strong Ba star HD 50082. The average masses are 1.9M and 1.5M respectively for mild Ba stars and strong Ba stars (Jorissen et al. 1998). We obtain masses of 1.90M and 1.32M for HD 27271 and HD 50082 respectively, which are consistent with the mass judgments well.

HD 26886 and HD 98839 also show over-abundant Ba and La, but the abundances are lower than those of HD 50082 and HD 27271. The elements Y, Zr of HD 26886 are weakly over-abundant. These results can be understood basically by the weaker overabundances of heavy elements in their AGB progenitor companions, by the lower neutron exposures occurred in the AGB progenitors.

## 6 Comparison with wind accretion model results

What is the reason to cause the heavy-element overabundances of Ba stars? It is generally believed that the overabundances are caused by binary accretion, namely, Ba stars accret the ejected materials from their companions, the former AGB stars (the present white dwarfs), which synthesized these heavy elements by themselves and ejected the elements into interstellar medium through stellar wind (Liang et al. 2000 and references therein).

In our previous paper (Liang et al. 2000), we calculated the s-process nucleosynthesis of AGB star with 3M and solar metallicity. At the same time, we set up an angular momentum conservation model of wind accretion in binary systems. Using such model we calculated the heavy-element abundances of Ba stars, and successfully explained the observed abundances of some samples taken from Za$`\stackrel{ˇ}{c}`$s (1994) and Busso et al. (1995).

We try to explain the abundance patterns of our sample stars using our theoretical method. The calculation of the heavy element abundances of Ba stars is made in two separate steps. Firstly, adopting the theory of s-process nucleosynthesis and the latest TP-AGB model (Straniero et al. 1995; Straniero et al. 1997; Gallino et al. 1998; Busso et al. 1999), we calculate the overabundances of the intrinsic AGB star at each ejection. Then, combining the accreting matter predicted by the model of wind accretion on successive occasions and mixing, we calculate the heavy-element overabundances of the barium star. The s-process nucleosynthesis scenario and the change equations of orbital semi-major axis $`A`$ and eccentricity $`e`$, $`\frac{\mathrm{\Delta }A}{A}`$ and $`\frac{e\mathrm{\Delta }e}{1e^2}`$, were given in Liang et al. (2000) and Liu et al. (2000).

The standard case of wind accretion is: $`M_{1,0}`$=3.0$`M_{}`$, $`M_{2,0}`$=1.3$`M_{}`$, $`v_{\mathrm{ej}}`$=15$`\mathrm{km}\mathrm{s}^1`$ ($`M_{1,0}`$ is the main sequence mass of the intrinsic AGB star, the present white dwarf, in the binary system; $`M_{2,0}`$ is the corresponding mass of the present Ba star; $`v_{\mathrm{ej}}`$ is the wind velocity). The used accretion rate is 0.15 times of the Bondi-Hoyle’s accretion rate (Liang et al. 2000; Boffin & Za$`\stackrel{ˇ}{c}`$s 1994).

We try to explain the observed abundances of our three Ba stars, HD 50082, HD 27271 and HD 26886, from theoretical calculations. Notice that there is no eccentricity of HD 98839 was observed, so we do not calculate the theoretical abundances to fit its observations. In our calculations, we try to make the calculated orbital period $`P`$ and eccentricity $`e`$ match the observations of the samples. Table 6 lists the observed orbital elements and their Ba classes of the samples taken from Jorissen et al. (1998), and their masses obtained by us (also in Table 1).

As the studies show, the Ba stars with orbital period $`P>`$1600 days can be formed through wind accretion (Zhang et al. 1999; Liang et al. 2000; Liu et al. 2000). Jorissen et al. (1998) suggested the corresponding period is 1500 days. Possibly, the Ba stars with lower orbital period form through other scenarios: dynamically stable late case C mass transfer or common envelope ejection.

For our sample of Ba stars, the orbital periods of HD 50082 and HD 27271 are 2896 and 1693.8 days ($`>`$1600 days) respectively. Their masses are consistent with the average masses of typical strong Ba stars (1.5M) and mild Ba stars (1.9M) (see Sect. 5.3), so we calculate the theoretical abundances of Ba stars using the standard case of our wind accretion model to compare with the observed patterns. Fig. 3(a), (b) show that the calculations (the solid lines) can fit observations very well within the errors for HD 50082 and HD 27271. The variable “$`a`$” represents the times of the corresponding standard exposures in the 13C profile suggested by Gallino et al. (1998) (details can be found in Liang et al. 2000). In summary, the higher $`a`$ value reflects the higher neutron exposure occurred in interiors of AGB progenitor companions.

Fig. 3(a), (b) show that the required neutron exposure characteristic parameter $`a`$ is higher for HD 50082, $`a`$=2.6; while for HD 27271, $`a`$=1.6 is adequate. This difference reflects the “strong” and “mild” Ba properties of each star.

For HD 26886, since its mass is 2.78M (though with large error), we adopt 2.5M as the initial mass of a Ba star for a binary system in our calculation, instead of 1.3M in the standard case. However, the fit failed (Fig. 3(c), the bottom panel). Maybe this is due to the lower orbital period of HD 26886 (1263.2 days) than the lower limit of wind accretion scenario (1600 or 1500 days). Perhaps dynamically stable late case C mass transfer or common envelope ejection scenario is the formation source of HD 26886. However, more studies are needed to understand this suggestion.

## 7 Discussions

### 7.1 Stellar masses and mild Ba stars, strong Ba stars

Adopting stellar evolution tracks given by Girardi et al. (2000) and some related parameters (see Sect. 3), we calculate the masses of the sample Ba stars (see Table. 1, 6). Their masses are consistent with the statistical and theoretical results. Table 9 of Jorissen et al. (1998) shows that the average mass of mild Ba stars is 1.9M with 0.60M companion white dwarfs. For HD 27271, a mild Ba star, its mass is 1.90M, which fits the average 1.9M of mild Ba stars very well. As for HD 26886, the calculated 2.78M is consistent with the average 2.3M of mild Ba stars with 0.67M companion white dwarfs (Jorissen et al. 1998) within the errors. HD 98839 is discussed in next section.

For the strong Ba star HD 50082, we obtain a mass of 1.32 M, which is consistent with the average mass 1.5M of strong Ba stars with 0.60M companion white dwarfs given by Jorissen et al. (1998).

Besides the different masses between the strong and mild Ba stars, their abundance patterns are different too (Sect. 5.3). What causes the discrepancy of abundances in these two Ba groups? Jorissen et al. (1998) carefully discuss the possible reasons. Considering our results, we emphasize and analyze two possible reasons again here: (i) a lower mass for the barium star results in a smaller dilution of the accreting materials in envelope of the Ba star; (ii) the lower metallicity results in larger heavy-element overabundances in the AGB progenitor. Our results are consistent with these two points. The strong Ba star HD 50082 has the lowest mass and lowest metallicity in the four samples. Thus it shows the highest heavy-element overabundances.

### 7.2 Inspiration from HD 98839: a “false Ba” star?

For HD 98839, the atmosphere parameters, \[Fe/H\] and over-rich Ba abundance obtained by us are similar to those of other researchers (Pilachowski 1977; Fernández-Villacañas et al. 1990; Mishenina 1996). We estimate its mass to be 3.62M. Its high mass and heavy-element overabundance pose some questions on Ba stars formation through binary accretion and also on nucleosynthesis scenario of AGB stars.

The mass of HD 98839 is 3.62M, higher than the average masses of mild Ba stars (1.9 or 2.3M, see Sect. 7.1). And there is also a discrepancy between this value and the suggestion of Han et al. (1995), in which the authors assigned masses of less than 3M for all Ba stars.

How can these results and discrepancies be understood? On one hand, this mild Ba star may have been classified erroneously as Ba star, namely it is “false Ba” star, and the real story is that it produced the overabundant heavy elements by itself. However, according to the derived parameters, $`M_{\mathrm{bol}}`$=$``$1.2 and $`T_{\mathrm{eff}}`$=4866 K, HD 98839 should not evolve to AGB stage. Generally, AGB stars have $`M_{\mathrm{bol}}`$ from $``$1.7 to $``$5.5 and $`T_{\mathrm{eff}}`$ from 3200 K to 3850 K respectively (Smith & Lambert 1985, 1986, 1990). In this case, if the heavy element overabundances of HD 98839 is an intrinsic nucleosynthesis event within itself, the over-rich s-process elements is usually associated with mixing process at the helium core flare, which is consistent with the alternative origin for the Ba ii star phenomenon suggested by Malaney (1987) and Malaney & Lambert (1988). Also, Mennessier et al. (1997) suggested that some of their mild Ba stars with mass higher than 3M may be lower luminosity, post He-flare stars.

Moreover, HD 98839 displays very long orbital period ($`>`$11000 days) up to the upper detected limit (Jorissen et al. 1998), and there is no observed binary orbital eccentricity. Thus, it is not unrealistic to consider it as an independent star.

However, if these high mass mild Ba stars are really “true Ba” stars, and the binary accretions cause their heavy element overabundances, it will challenge the present nucleosynthesis scenario of AGB stars. At present, the popular viewpoint is that 1$``$3M low mass AGB stars with 13C neutron source are the main nucleosynthesis sites of s-process elements, while the intermediate mass AGB stars with 22Ne neutron source cannot fit to the observations (Malaney 1987; Busso et al. 1995; Busso et al. 1999; Liang et al. 2000). However, if masses of Ba stars are higher than 3M, their companions should be intermediate mass stars on the main-sequence. How the intermediate mass AGB stars result in the heavy element overabundance of these Ba stars? If this is the case, what kind of neutron exposure scenario works? Single exposure? More studies are needed on s-process nucleosynthesis scenario of AGB stars. Certainly, we need also the large sample of high mass mild Ba stars to check its abundance patterns.

In addition to the high mass (3.62M) of mild Ba star, HD 98839, obtained by us, there are also some other high mass samples. Smith et al. (1980) calculated mass of $`\zeta `$ Cap is about 3.5M, but with large error (+3.6,–1.8M) (also see Boffin & Jorissen 1988). Mennessier et al. (1997) studied higher mass Ba stars. They estimated that the masses are up to 7M for their group S and 4.5M for group C in which there are mild Ba stars. Monte Carlo simulations of Karakas et al. (2000) predict that the masses of Ba stars may be up to 6M.

Therefore, it is necessary to pay more attention on these high mass mild Ba stars.

## 8 Summary

We obtain abundances of $`\alpha `$, iron peak and s-process elements in four Ba stars, HD 50082, HD 27271, HD 26886 and HD 98839, using high resolution and high S/N spectra. The results show that all of samples are disk stars. They show the similar metallicities and the similar abundances of $`\alpha `$-elments and iron peak elements to the solar values. All of them show s-process element overabundances relative the solar. We obtain their masses using the highly precise Hipparcos data, precise photometric parameters and stellar evolution tracks. The mass results fit well the average characteristic masses of strong and mild Ba stars. And the high mass of mild Ba star HD 98839 pose more thoughts on Ba star phenomenon.

Using our wind accretion model of binary systems and AGB stars s-element nucleosynthesis calculation, we conclude that the wind accretion formation scenario is possible for HD 50082 and HD 27271. Wind accretion cannot explain the observations of HD 26886 with the short orbital period ($`P<`$1600 days). The mild Ba star HD 98839 with high mass and very long orbital period may be either a “false Ba” star with the heavy elements enriched by itself or a “true Ba” star.

Fortunately, though we only have four samples, they represent wide Ba star groups. From stellar masses, there is one strong Ba star among them, and also mild Ba stars; there are the members with typical mass of Ba stars, also there are the high mass mild Ba stars, which are controversial. From orbital elements, there are samples in typical wind accretion orbital period range, there is one sample with very long orbital period, up to the detected limit, and there is one sample with shorter value, which may correspond to other formation scenarios. So these samples are very interesting in understanding the abundance patterns and formation scenarios of Ba stars.

Certainly, we need larger samples to understand the nature of Ba stars. At present, there is a large sample of Ba stars with orbital elements (Carquillat et al. 1998; Udry et al. 1998a, 1998b; Jorissen et al. 1998), but without the corresponding element abundances. It is very interesting and necessary to obtain their element abundances to be combined with the orbital elements. It is specially important to compare the related results of the two groups with different orbital periods: $`P<`$1600 and $`P>`$1600 days, and to study the stars with very long orbital periods. These results will help us to understand the various formation scenarios of Ba stars. At the same time, it is important to obtain their masses, particularly, for the possible high mass mild Ba stars. Maybe these high mass mild Ba stars can give more clues on AGB stars nucleosynthesis. Moreover, it will be helpful to understand the different abundance patterns between “strong” and “mild” Ba stars.

There is a larger sample of Ba stars (318 from L$`\ddot{u}`$ catalogue 1991) have been studied whose absolute magnitudes and kinematics by Gómez et al. (1997). Their results show that the samples of Ba stars are an inhomogeneous group, three of them among the four groups are disk population stars, and the fourth group contains halo stars. Mennessier et al. (1997) confirmed further this inhomogenity. If we can supplement the heavy-element abundances to compare the various characters of these different groups, it will be a stronger method to understand the natural properties of Ba stars.

In addition, it is better to obtain the carbon and oxygen abundances together with s-process elements. The reason is that AGB stars are important nucleosynthesis sources of element carbon, which is taken out from stellar interiors together with s-process elements by the third dredge-up process (Liang et al. 2000; Liang et al. 2001). Element carbon will be then accreted by Ba stars, the companions of binary systems.

Considering Ba stars from stellar population, it is interesting to combine Ba stars with CH stars. Some reseachers suggest that metal-rich stars evolve into Ba stars and the metal-poor stars evolve to CH stars, which are metal-poor stars (Vanture 1992a, 1992b). At present, it is believable that nucleosynthesis of metal-poor AGB stars are more efficient than the metal-rich one (Busso et al. 1999; Busso et al. 2001). Thus, it is valuable to study the heavy-element abundances of metal-poor CH stars. But CH stars are fainter, some are as low as 12 mag, which makes more difficult to get the required spectra.

CH subgiants (Luck & Bond 1982, 1991; Smith et al. 1993), Ba dwarfs (North et al. 1994) and extrinsic S stars (Jorissen et al. 1998) should be the closest groups to Ba giants. It will be interesting to combine and compare all of these group stars. Perhaps we will then be able to study Ba stars in an evolutionary sequence.

###### Acknowledgements.

We thank the anonymous referee for very useful suggestions on the original manuscript. We thank Dr. Rafael Guzmán for his great help for improving the English expression. Y. C. Liang thanks Dr. Francois Hammer and other staffs of GEPI for their friendly help. Y. C. Liang thanks Huawei Zhang, Jianrong Shi, Zhenxi Zhu, Guangshun Bai and Gang Liu for their help with data reduction. This research work is supported by the National Natural Science Foundation of China under grant No. 10173014 and NKBRSF G1999075406.


text:
# Selection of Tau Leptons with the CDF II Trigger System

## 1 Introduction

Run II of the Tevatron produces $`p\overline{p}`$ interactions with an instantaneous luminosity about $`10^{32}`$ $`\mathrm{cm}^2\mathrm{s}^1`$. The CDF Detector trigger system  decreases the resulting initial event rate of 2.5 MHz down to 70 Hz for data storage and offline analysis, with the goal of efficient extraction of interesting physics events from the large minimum bias background. One important subset of these physics signatures includes $`\tau `$ leptons. Events with $`\tau `$’s are important for the study of various Standard Model (SM) processes such as $`Z\tau \tau `$, $`t\overline{t}`$ with $`W\tau \nu _\tau `$, $`H\tau \tau `$ and also to extend the reach of searches for physics beyond the SM. These physics goals drive us to look for $`\tau `$-like objects starting with Level-3 of the trigger, and our implementation selects events with two leptons in the final state: $`e\tau _h,\mu \tau _h,\tau _h\tau _h`$ and also $`\tau _h\nu `$, as well as $`ee,e\mu `$ and $`\mu \mu `$. With ”$`\tau _h`$” we indicate hadronically decaying $`\tau `$’s, while the $`e`$ or $`\mu `$ can be either produced directly or through a leptonic decay of the $`\tau `$.

These tau triggers benefit from the upgraded CDF trigger system, in particular the tracking processor at Level-1 (L1) and the refined tracking information available at Level-2 (L2) .

## 2 Implementation of the $`\tau `$ Triggers

### 2.1 $`\tau `$ reconstruction

Taus promptly decay in leptonic or hadronic (65%) modes with at least one $`\nu `$ in the final state. Hadronic $`\tau `$ decays have the distinct signature of a narrow isolated jet with low multiplicity (1 or 3 prongs) and low visible mass ($`<M_\tau `$). The basis for our triggers, therefore, is a $`\tau `$-cone algorithm for the reconstruction of these decays. To build the $`\tau `$-cone object, we start with a narrow calorimeter cluster, above a suitable $`E_T`$ threshold, matched to a seed track with momentum above a $`P_T`$ cut. The region within an angle $`\mathrm{\Theta }_{Sig}`$ from the seed track direction is used to define a cone of tracks to be associated with the $`\tau `$. The region between $`\mathrm{\Theta }_{Sig}`$ and $`\mathrm{\Theta }_{Iso}`$ defines the isolation cone: we require that no tracks with $`P_T`$ higher than a fixed low threshold be found in this region. At the Level-3 (L3) of the trigger, $`\mathrm{\Theta }_{Sig}`$ and $`\mathrm{\Theta }_{Iso}`$ have values of $`10^{}`$, and $`30^{}`$, respectively, and are determined in 3 dimensions.

We have implemented 4 different triggers for physics processes with $`\tau `$’s in the final state, installed in the CDF trigger tables in January 2002. Below we describe their main features.

### 2.2 Electron plus track, Muon plus track

The selection for the $`e`$ plus track trigger starts with L1 and L2 requirements of a single EM tower with $`E_T>`$ 8 GeV with an associated XFT track with $`P_T>`$ 8 GeV/$`c`$ and a second track with $`P_T>`$ 5 GeV/$`c`$. At L3, these conditions are refined and charged track isolation around the reconstructed (2nd) track is imposed. The current average cross-section for this L3 trigger is $``$29 nb. For the muon plus track trigger we require a muon stub matched at L1 to an XFT track with $`P_T>\mathrm{\hspace{0.33em}4}GeV/c`$, with an increased track threshold of 8 GeV/$`c`$ at L2. The ”track” requirements are identical to those for the electron plus track trigger. The current average cross section for this trigger is $``$16 nb.

### 2.3 Di-$`\tau `$,$`\tau `$ plus Missing Transverse Energy

At L1 the di-$`\tau `$ trigger requires 2 calorimeter towers with $`E_T>`$ 5 GeV and 2 matching XFT tracks with $`P_T>`$ 6 GeV/$`c`$, separated by an angle of $`\varphi >30^{}`$. L2 requires cluster $`E_T>`$ 10 GeV and imposes track isolation. At L3, using the full reconstruction code, 2 $`\tau `$ candidates with seed track $`P_T>`$ 6 GeV/$`c`$, and originating from the same vertex: $`|\mathrm{\Delta }Z|<10`$ cm, are required. The L3 cross section for this trigger is $``$12 nb.

The $`\tau `$ plus missing transverse energy (MET) trigger requires $`E_{MET}>`$10 GeV at L1. At L2 this is increased to $`E_{MET}>`$20 GeV and a calorimeter cluster with an isolated track with $`P_T>`$10 GeV/$`c`$ is required. This is followed by full event reconstruction at L3, requiring at least one $`\tau `$ candidate having seed track $`P_T>`$ 4.5 GeV/$`c`$. The present trigger cross section is 5 nb.

## 3 First Results

The study of $`Z\tau _e\tau _h`$ is performed using $`72\mathrm{pb}^1`$ of data taken with the electron plus track trigger . The track multiplicity associated with the $`\tau _h`$ candidate found in events passing selection cuts is shown in Fig. 1. We observe a clear $`\tau `$ signal over background levels even before requiring an opposite sign charge for the electron and the $`\tau _h`$. The mass distribution of the electron-$`\tau _h`$-MET system in the data is also found to be consistent with the $`Z\tau \tau `$ hypothesis.


text:
# Approach to Vibrational to Axially Rotational Shape Phase Transition Along the Yrast Line

## Abstract

With the energy surface of the nucleus in U(5) symmetry being analyzed in the framework of thermodynamics, the vibration and rotation phase diagram in terms of the angular momentum and deformation parameter is given. Together with examining the energy spectrum, we propose a theoretical approach to describe the vibrational to axially rotational phase transition along the yrast line. By analyzing the available experimental data we show that the vibrational to rotational shape phase transition along the yrast line takes place in many nuclei.

PACS No. 21.10.Re, 21.60.Fw, 27.50.+e, 27.60.+j, 27.70.+q

It has been well known that shape phase transition is one of the most significant topics in nuclear structure research. Many evidences of nuclear shape phase transition have been observed. For instance, in several isotopes, there exists shape phase transition from vibration to axial rotation or $`\gamma `$-unstable rotation with respect to the variation of neutron number, and a triple point may appear. Even in one mode of collective motion there may involve different characteristics, for example, along the yrast line there exists transition between rotations holding different relations between angular momentum and rotational frequency (referred to as band crossing and exhibiting a backbending). Very recently it was found that there involves evolution from vibration to rotation along the yrast line of a nucleus. On the theoretical side, the interacting boson model (IBM) has been shown to be successful in describing the shape phase transition along a series of isotopes. And analytic solutions for the critical points of the phase transitions have been found. The cranked shell model (CSM) has been known to be able to describe the band crossing very well. However, a theoretical approach to describe the shape phase transition from vibration to rotation along the yrast line in a nucleus has not yet been established. Especially the phase diagram in terms of the angular momentum and the deformation parameter has not yet been obtained. Even so, the energy level statistics analysis has shown that the states in the U(5) symmetry with $`A<0`$ involve the coexistence of vibrational and axially rotational shapes. By analyzing the energy surface and the energy spectrum of the U(5) symmetry, we will show that the U(5) symmetry with a special choice of parameters can be an approach to describe the shape phase transition along the yrast line and such kind phase transition takes place in many nuclei.

In the original version of the IBM (IBM1), the collective states of nuclei are described by the coherent states of $`s`$\- and $`d`$-bosons. The corresponding dynamical group is U(6), and involves three dynamical symmetry limits U(5), O(6) and SU(3). Taking into account one- and two-body interactions among the bosons, one has the Hamiltonian of the nucleus in U(5) symmetry as

$$\begin{array}{c}\text{ }H_{U(5)}=E_0+\epsilon _dC_{1U(5)}+AC_{2U(5)}+BC_{2O(5)}+CC_{2O(3)},(1)\text{ }\hfill \end{array}$$

where $`C_{kG}`$ is the $`k`$-rank Casimir operator of the group $`G`$ and the parameters hold relation $`C|B||A|\epsilon _d`$. Making use of the coherent state formalism of the IBM, one can express the energy functional for the U(5) Hamiltonian as

$$\begin{array}{c}\text{ }E(N,\epsilon ,A;\beta )=E_0+\epsilon N\frac{\beta ^2}{1+\beta ^2}+AN(N1)\frac{\beta ^4}{(1+\beta ^2)^2},(2)\text{ }\hfill \end{array}$$

where $`\beta `$ is a monotonous function of the usual deformation parameter (for example, for the rare earth nuclei, the usual quadrupole deformation parameter $`\beta _2`$ and the $`\beta `$ hold relation $`\beta _20.16\beta `$) and the $`\epsilon `$ can be expressed in terms of the parameters in Eq. (1) as $`\epsilon =\epsilon _d+5A+4B+6C`$. If $`\beta ^2`$ is small (less than one), Eq. (2) can be rewritten as

$$\begin{array}{c}\text{ }E(N,\epsilon ,A;\beta )=E_0+\epsilon N\beta ^2+[AN(N1)\epsilon N]\beta ^4+[\epsilon N2AN(N1)]\beta ^6+\mathrm{}\mathrm{},(3)\text{ }\hfill \end{array}$$

On the basis that the “deformation parameter” $`\beta `$ is a quantity to characterize the shape of a nucleus, one can take the $`\beta `$ as the order parameter to identify the nuclear shape phase transition.

In the general framework of Landau’s phase transition theory, it has been pointed out that, if the Landau free energy $`F(T,\xi )`$ can be expressed as

$$\begin{array}{c}\text{ }F(T,\xi )=F_0+\frac{1}{2}\alpha (TT_0)\xi ^2+\frac{1}{4}g_4\xi ^4+\frac{1}{6}g_6\xi ^6+\mathrm{}\mathrm{}(4)\text{ }\hfill \end{array}$$

where $`T`$ is the control quantity, $`\xi `$ is the order parameter, $`\alpha `$ is a constant. It has been known that, if the other coefficients which depend also on the control quantity $`T`$ hold relation $`g_4<0`$, $`g_6>0`$, the system may involve a first order phase transition.

For a nucleus involving vibration to rotation phase transition along the yrast line, the control quantity is the angular momentum $`L`$. Comparing Eq. (3) with Eq. (4), we propose that the Eq. (3) can be rewritten as

$$\begin{array}{c}\text{ }E(N,\epsilon ,A;\beta )=E_0+\frac{1}{2}\alpha (LL_0)\beta ^2[\frac{1}{2}\alpha (LL_0)AN(N1)]\beta ^4\text{ }\hfill \\ \text{ }+[\frac{1}{2}\alpha (LL_0)2AN(N1)]\beta ^6+\mathrm{}\mathrm{},(5)\text{ }\hfill \end{array}$$

where $`AN(N1)`$ relies also on the control quantity $`L`$. It is evident that, if $`\alpha <0`$ (i.e., $`\alpha (LL_0)>0`$ for $`L<L_0`$) and $`A<0`$, Eq. (5) stands for the free energy of a nucleus just the same as Eq. (4) with $`g_4<0`$, $`g_6>0`$. Even if $`A>0`$, we may have $`\frac{1}{2}\alpha (LL_0)2AN(N1)>0`$, Eq. (5) corresponds also to Eq. (4) with $`g_4<0`$ and $`g_6>0`$. Analyzing the free energy in Eq. (5) more cautiously, one can recognize that there exists an angular momentum $`L_0`$ with which the $`\frac{^2E(N,\epsilon ,A;\beta )}{\beta ^2}=0`$. And $`\alpha (LL_0)>0`$, $`\frac{^2E(N,\epsilon ,A;\beta )}{\beta ^2}>0`$, if $`L<L_0`$; $`\alpha (LL_0)<0`$, $`\frac{^2E(N,\epsilon ,A;\beta )}{\beta ^2}<0`$, if $`L>L_0`$. Meanwhile, there exists an angular momentum $`L_c=L_0(\frac{4}{3}\sqrt{3}2)\frac{AN(N1)}{\alpha }`$ (if $`A<0`$) with which the free energy involves one maximum and two equal minima, one of which is generated at $`\beta =0`$, another corresponds to $`\beta 0`$. There exists also a maximal angular momentum $`L_{max}=L_0+(\sqrt{6}+2)\frac{AN(N1)}{\alpha }`$ (if $`A<0`$). If $`LL_{max}`$, the energy surface has only one maximum at $`\beta =0`$ but no minimum. The $`\beta `$ dependence of the free energy at some typical angular momenta is illustrated in Fig. 1. The figure shows apparently that, as the angular momentum $`L<L_c`$ the nucleus has a stable vibrational phase (the free energy takes the smaller minimum at $`\beta =0`$). If the angular momentum $`L_cL<L_0`$, the nucleus has a stable rotational phase (the free energy takes the smaller minimum at $`\beta 0`$). Meanwhile there may involve coexistence of rotation and vibration, since the energy surface holds two minima at $`\beta =0`$, $`\beta 0`$, respectively. If the angular momentum $`L_0<L<L_{max}`$, the nucleus has only a rotational phase. The angular momentum $`L_c`$ is definitely the critical angular momentum for the phase transition from vibration to rotation to happen. In addition, if $`LL_{max}`$, the nucleus will collapse. We obtain then explicitly how the stable phase of a nucleus is determined by the angular momentum and the “deformation parameter”. The phase diagram in terms of the angular momentum $`L`$ and the “deformation parameter” $`\beta `$ can thus be displayed in Fig. 2. We reach thus a conclusion that the U(5) symmetry with $`A<0`$ can describe well the vibrational to rotational nuclear shape phase transition along the yrast line. And the transition is a first order phase transition. Even if $`A>0`$, the first order phase transition may also happen in the system if the interaction strengths satisfy the constraint mentioned above. It manifests that, the vibrational to rotational phase transition along the yrast line can take place in the nucleus holding the U(5) symmetry in IBM1.

On the phenomenological side, the energy spectrum of a nucleus in U(5) symmetry can be given as

$$\begin{array}{c}\text{ }E_{U(5)}=E_0+\epsilon _dn_d+An_d(n_d+4)+B\tau (\tau +3)+CL(L+1).(6)\text{ }\hfill \end{array}$$

where $`n_d`$, $`\tau `$, and $`L`$ are the irreducible representations (IRREPs) of the group U(5), O(5) and O(3), respectively.

From the spectrum generating process one can recognize that, if $`A>0`$ and $`B>0`$, the states with $`\tau =n_d`$, $`L=2n_d`$ form the ground state band and the yrast band, and appear as the anharmonic vibrational ones with increasing frequency $`\mathrm{}\omega =\epsilon _d+(n_d+4)A`$. If $`A<0`$, the $`E_{gsb}(n_d)`$ ($`E_{gsb}(L)`$) is an upper-convex parabola against the $`n_d`$ ($`L`$), and can describe the collective backbending of high spin states. In the case with $`A<0`$, $`B<0`$, the freedom of $`\tau `$ in Eq. (6) can be washed out. It turns out then the state with $`\tau =n_d`$ may still be the yrast one and there exists a d-boson number $`n_d^{(c)}`$ and an angular momentum

$$\begin{array}{c}\text{ }L_c=2n_d^{(c)}=\frac{2(\epsilon _d+4A+3B)}{A+B}2N_0,(7)\text{ }\hfill \end{array}$$

where $`N_0=N`$ with $`N`$ being the total boson number. As the angular momentum $`LL_c`$, the yrast states are no longer the anharmonic vibrational ones mentioned above, but the quasi-rotational ones with $`\tau =n_d=N_0`$. It indicates that the U(5) symmetry with parameters $`A<0`$ and $`B<0`$ may describe the vibrational to rotational phase transition along the yrast line. The $`L_c`$ given in Eq. (7) is the critical angular momentum. For each yrast state with angular momentum $`LL_c`$, its energy can be given as

$$\begin{array}{c}\text{ }E(LL_c)=E_0+\frac{A+B+4C}{4}L^2+\frac{\epsilon _d+4A+3B+2C}{2}L.(8)\text{ }\hfill \end{array}$$

and can be illustrated as a part of an upper-convex parabola. Whereas for the one with angular momentum $`L>L_c`$, its energy should be expressed as

$$\begin{array}{c}\text{ }E(L>L_c)=E_0^{}+CL(L+1),(9)\text{ }\hfill \end{array}$$

and can be displayed as a part of an upper-concave parabola. For instance, for a system with $`N=15`$ and parameters $`\epsilon _d=0.80`$ MeV, $`A=0.025`$ MeV, $`B=0.01`$ MeV and $`C=0.004`$ MeV, with Eq. (7) we can fix the critical angular momentum $`L_c=8\mathrm{}`$. The energy of the yrast states against the angular momentum $`L`$ can be illustrated in the left panel of Fig. 3. It has been known that the energy of E2 transition $`\gamma `$-ray over spin (E-GOS) $`R=\frac{E_\gamma (LL2)}{L}`$ can be taken as a quite good signature to manifest the vibrational to axially rotational phase transition along the yrast line. As an auxiliary evidence, we show also the E-GOS of the yrast states with these parameters in the right panel of Fig. 3. The figure indicates apparently that the yrast states involve a vibrational to axially rotational phase transition and the angular momentum $`L_c`$ is definitely the critical point for the phase transition to take place. Such a transition is quite similar to that between the states with $`n_p=0`$ and $`n_p=N`$ in the vibron model with random interactions.

In the case $`A<0`$, $`B>0`$ or $`A>0`$, $`B<0`$, there may exist a critical angular momentum $`L_c`$ with $`\frac{L_c}{2}N_0<N`$, too. For the states with $`LL_c`$, if $`\epsilon _d+A(N_0+4)+B(N_0+3)>0`$, their energies behave like an anharmonic vibrator’s. For the ones with $`L>L_c`$, their energies behave as the quasi-rotational states with $`N_0`$ bosons. It is also possible that the $`N_0`$ takes several values less than the $`N`$. Then the second and more higher backbending may emerge.

We have analyzed the available experimental data of the yrast states of the nuclei with $`30Z100`$ and simulated the data with least-square fitting in our present approach. Both the experimental data and the best fitted energy spectra of the yrast states in some even-even nuclei are illustrated in Fig. 4 (the other data are available if required). It shows obviously that, besides the ones identified in the $`A110`$ mass region, there exist vibrational to axially rotational phase transition along the yrast line in other nuclei. For example, 72Se, 90Zr, 92Zr, 98Pd, 108Sn, 116Sn, 136Te, 146Gd, 148Gd, 148Dy, 150Dy, 186Hg, 194Pb, 196Pb, 198Pb, 200Po, 202Po, 206Po, 214Rn and others are the ones whose anharmonic vibrational frequency in the corresponding phase decreases against the increasing of angular momentum. While 66Zn, 66Ge, 68Se, 82Sr, 84Sr, 86Mo, 114Pd, 116Pd, 116Cd, 118Cd, 120Cd, 112Te, 114Te, 116Te, 116Xe, 118Xe, 124Xe, 128Xe, 130Xe, 124Ba, 130Ba, 148Sm, 142Gd, 144Dy, 156Er, 158Yb, 160Yb, 160Hf, 162Hf, 164W, 170Os, 188Pt, 192Pt, 194Pt, 184Hg and so on are the ones whose vibrational frequency increases against the increasing of angular momentum.

In summary, by analyzing the energy surface of the nucleus in U(5) symmetry in the framework of thermodynamics, we give a phase diagram of the vibration and the rotation in terms of the angular momentum and the deformation parameter in this letter. We have then proposed an approach to describe the vibrational to axially rotational phase transition along the yrast line. With the energy spectrum being examined, we show that our presently proposed approach can describe very well the vibrational to axially rotational phase transition along the yrast line. We have also analyzed the energy spectra of the even-even nuclei with $`30Z100`$. It shows that, besides the ones in $`A110`$ mass region identified by Regan and collaborators, there exist many other nuclei involving vibrational to axially rotational phase transition along the yrast line.

This work is supported by the National Natural Science Foundation of China under the contract Nos. 19875001, 10075002, and 10135030 and the Major State Basic Research Development Program under contract No. G2000077400. One of the authors (Y.X. Liu) thanks the support by the Foundation for University Key Teacher by the Ministry of Education, China, too.


text:
# Structural and magnetic properties of transition metal substituted ZnO

## I Introduction

A widegap II-VI semiconductor, ZnO attracts attention as a material with possible application in optoelectronic devices such as solar cells, ultraviolet emitting diodes and transparent high-power electronic devices. The discovery of ferromagnetism at temperatures above 100 K in the III-V semiconductor Ga1-xMnxAs Ohno98  provided practical means to incorporate spin into semiconductor electronics. Theoretical predictions of room temperature ferromagnetism in diluted magnetic semiconductors Dietl00  recently focused attention on magnetic-ion-substituted ZnO with a wurtzite structure similar to GaAs. According to these calculations, room temperature ferromagnetism can exist in $`p`$-type doped Zn1-xMnxO with 5%Mn and $`3.510^{20}`$ holes per cm3. Ab initio band calculations Sato00a  predict stability of ferromagnetism in $`p`$-type Zn1-xMnxO, and antiferromagnetism in $`n`$-type Zn1-xMnxO. Similar calculations predict a ferromagnetic phase for both carrier-undoped and $`n`$-type ZnO substituted with Fe, Co, or Ni. Sato00b  The parent compound ZnO without intentional carrier doping shows $`n`$-type conduction related to oxygen vacancies and Zn interstitials. Minegishi97  The introduction of $`p`$-type conduction in ZnO is difficult. Several methods, including nitrogen doping Minegishi97 ; Look02  and Ga and N codoping Yamamoto99 ; Joseph99 , have been reported. Successful preparation of $`p`$-type ZnO could allow fabrication of transparent $`pn`$ junctions and could enable synthesis of substituted ferromagnets.

Pulsed-laser deposited Zn1-xMnxO thin films with up to 35% Mn show spin-glass behavior. Fukumura01  We have previously shown for polycrystalline Zn1-xMnxO that the spin-glass behavior can be induced by high-pressure oxygen annealing and it is related to precipitation of impurity phase ZnMnO3. Kolesnik02  ZnO films doped with Co have been reported to be ferromagnetic with Curie temperature of about 280 K. Ueda01  However, Kim et al. pointed to the presence of Co clusters in their thin films as a reason for room temperature ferromagnetism. Kim02  High-temperature ferromagnetism in thin films of Co-doped anatase TiO2 Matsumoto01  has also been explained by the fact that cobalt is not substitutional for titanium in TiO2. Stampe02  In this study, we investigate polycrystalline Zn1-xTMxO, where TM = Mn, Fe, and Co. We show that no high-temperature ferromagnetic ordering is present in single-phase samples. Ferromagnetism in these materials can be induced by second phases that may appear after non-optimal synthesis conditions.

## II Experimental methods

The Zn1-xTMxO samples in this study were prepared using a standard solid-state reaction technique. Mixtures of ZnO and MnO2, Fe2O3, or Co3O4, for TM = Mn, Fe, and Co, respectively, were fired in air several times at elevated temperatures 900, 1000, 1100, and 1200C with intermediate grindings. The samples were then annealed in an atmosphere of various gases (Ar, H2, 1%H2/Ar, and O2 under high pressure of 600 bar). High pressure oxygen annealing was intended to produce $`p`$-type doped ZnO materials. Magnetic ac susceptibility and dc magnetization were measured using a Physical Property Measurement System and a Magnetic Property Measurement System (both Quantum Design) at temperatures up to 400 and 800 K, respectively. Energy dispersive x-ray spectroscopy (EDXS) analysis was performed by a Hitachi S-4700-II scanning electron microscope. Typically, 15-20 EDXS spectra were collected for each composition at various locations across the surface of a sintered pellet. Using the average Zn and TM cation contents, which were obtained from the EDXS spectra we determined the effective composition Zn$`_{1x_{\mathrm{eff}}}`$TM$`_{x_{\mathrm{eff}}}`$O of the studied samples. X-ray diffraction experiments have been performed using a Rigaku x-ray diffractometer. We have analyzed the x-ray diffraction patterns using the Rietveld technique with the General Structure GSAS code.GSAS85  The refinements were done using the wurtzite hexagonal space group P63mc with (Zn,TM) atoms located at (2/3, 1/3, 0) and O atoms located at (2/3, 1/3, $`u`$). For both single-phase and multi-phase samples (where impurity peaks were observed), this was the only phase refined.

## III Relationship between synthesis conditions and structural properties

Single-phase samples can be synthesized in air for TM = Mn, and Co and in 1%H2/Ar for TM = Fe. X-ray diffraction data show that the Zn1-xTMxO samples are single-phase with the wurtzite structure up to a certain value of $`x`$. We estimated this value (the solubility limit) to be equal to 0.1, 0.1, and 0.2 for TM = Mn, Fe, and Co, respectively. Above the solubility limit, characteristic impurities can be observed by x-ray diffraction, namely ZnMn2O4, FezO or (Zn,Fe)3O4, and CoO for TM = Mn, Fe, and Co, respectively.

In Fig. 1 X-ray diffraction patterns are presented for Zn1-xFexO after synthesis in 1%H2/Ar at 1000C. The samples with $`x0.1`$ are single-phase, and the presence of the second phase FezO can be observed for $`x0.15`$. Wüstite FezO can typically be obtained in a slightly non-stoichiometric form with $`z<1`$ depending on synthesis conditions; however, due to a very small amount of this phase in the Zn0.85Fe0.15O sample, we could not determine the value of $`z`$. Thus, we have separately prepared wüstite FezO by reduction of Fe2O3 in the same atmosphere, 1%H2/Ar at 1050C, and obtained the compound with the rock-salt structure and lattice constant $`a=4.3101`$ Å that is consistent with $`z=0.95`$. Stability of Zn1-xFexO and second phases are sensitive to changes of the synthesis temperature. The annealing of these samples in 1%H2/Ar at 950C leads to precipitation of second phase (Zn,Fe)3O4, which accounts for room temperature ferromagnetism in multi-phase Zn1-xFexO. Substitution of Fe in ZnO is not possible in air or at high pressure of oxygen at any temperature.

Fig. 2 shows x-ray diffraction patterns for the Zn0.8Co0.2O sample obtained from several annealing conditions. Various impurity peaks appear after annealing under reducing or strongly oxidizing conditions. These peaks are marked with stars in Fig. 2. After argon annealing at $`T=1100^{}`$C, rock-salt-like antiferromagnetic CoO impurity is present. A slight reduction of Zn0.8Co0.2O in 1%H2/Ar at $`1000^{}`$C leads to precipitation of Co metal, which causes a ferromagnetic contribution to the measured ac susceptibility and dc magnetization. High-pressure oxygen annealing at $`T=1080^{}`$C induces spinel-like ZnCo2O4 impurity. High-pressure oxygen annealing was found, thus, to favor formation of compounds other than Zn1-xTMxO, for which TM ions display oxidation states higher than 2+.

Structural parameters of Zn1-xTMxO obtained from the Rietveld refinements are presented in Fig. 3. Open symbols in Fig. 3 denote the multi-phase compositions, i.e. the ones for which the presence of impurity peaks is observed in the X-ray diffraction data. The unit cell volume changes almost linearly with substitution of TM for single-phase compositions. This behavior is expected according to the difference between ionic radii $`r`$ of tetrahedrally coordinated Zn2+ $`r=0.60`$ Å and other transition metal TM2+ ions ($`r`$ = 0.66, 0.63, and 0.58 Å for TM = Mn, Fe, and Co, respectively). Shannon76  The observed slight increase of the cell volume in case of TM = Co is inconsistent with the smaller ionic radius of Co than Zn estimated by Shannon. Shannon76  This increase can not be explained either by the presence of Co3+ in Zn1-xCoxO as Co3+ is smaller than Co2+. Shannon76  Apparently, the ionic size of tetrahedrally coordinated Co2+ ion is larger than that of Zn2+. The $`u`$ parameter ($`z`$-coordinate of the oxygen atoms) of the wurtzite structure of Zn1-xTMxO initially increases with increasing $`x`$, reaches its maximum close to the solubility limit, and subsequently decreases. The $`c/a`$ ratio \[Fig. 3(d)\] decreases slightly at a similar rate for TM = Mn and Co, while it decreases much more rapidly for TM = Fe. This effect shows that Fe-doped ZnO is more anisotropic than ZnO doped with other transition metal ions. Ref. Schulz79  has shown a phenomenological linear dependence $`u(a/c)^2`$ for all undoped wurtzite-like compounds including ZnO. Our observation for transition metal doped ZnO is qualitatively consistent with this tendency up to the solubility limit. The coefficient of this linear dependence (varying from 9 to 34) is much higher than the reported value of 1/3. Schulz79

The effective TM content $`x_{\mathrm{eff}}`$ in Zn1-xTMxO determined from EDXS analysis is presented in Fig. 4. In our previous work Kolesnik02  we have shown (based on the assumption that the ac susceptibility reflects spin-only behavior) that $`x_{\mathrm{eff}}<x`$ for TM = Mn (solid line in Fig. 4.) Present EDXS data for TM = Mn confirm that result. For TM = Co, we observe a good agreement between the effective and nominal TM contents, within the experimental error. For TM = Fe, $`x_{\mathrm{eff}}x`$ for $`x=0.05`$. However, for higher Fe concentration, $`x=0.1`$, we observe $`x_{\mathrm{eff}}0.14`$. This larger effective Fe content than the starting content, $`x`$, could arise from Zn volatilization during synthesis under reducing conditions in 1%H2/Ar. We did not include this composition for the analysis of the magnetic properties of Zn1-xFexO.

## IV Magnetic properties

Magnetic susceptibility is presented in Fig. 5 for Zn0.8Co0.2O. Here and throughout this paper, the diamagnetic ac susceptibility of $`0.3310^6`$ emu/g for ZnO Handbook00  was subtracted from the measured magnetic susceptibility. The ac susceptibility of the TM-doped ZnO resembles Curie-Weiss behavior which is also characteristic for other TM-containing semimagnetic semiconductors. Spalek86  The annealing in Ar under atmospheric pressure (which is correlated to the appearance of antiferromagnetic CoO impurity \[Fig. 2(b)\]) does not significantly change the value of the ac susceptibility. A significant decrease of the ac susceptibility can be observed after high-pressure oxygen annealing. After this annealing, a spinel-like ZnCo2O4 second phase \[Fig. 2(c)\] is present in the sample. The most dramatic change of the magnetic properties of Zn1-xCoxO can be observed after annealing in 1%H2/Ar at 1000C. The sample becomes ferromagnetic at a temperature higher than 800 K. This ferromagnetic behavior is a result of the presence of Co-metal impurity in the compound \[Fig. 2(d)\].

Inverse magnetic susceptibility for Zn1-xCoxO is shown in Fig. 6. The behavior presented here is characteristic to all studied TM-doped ZnO samples. These results, for TM = Mn, were presented elsewhere. Kolesnik02  We have analyzed ac susceptibility results within the framework of the diluted Heisenberg antiferromagnet theory of Spałek et al. Spalek86  Although this theory was developed for larger group VI elements, Se, Te, and S, we use it here also for the oxygen ion. At higher temperatures, the ac mass susceptibility can be described by the formula

$$\chi =\frac{C_M(x)}{T\mathrm{\Theta }(x)},$$
(1)
where $`\mathrm{\Theta }(x)=\mathrm{\Theta }_0x`$ is the Curie-Weiss temperature, $`C_M(x)=C_0x`$ is the molar Curie constant and $`C_0`$ is defined as

$$C_0=\frac{N(g_{\mathrm{eff}}\mu _B)^2S(S+1)}{3k_B\rho },$$
(2)
$`N`$ is the number of cations per unit volume, $`g_{\mathrm{eff}}`$ is the effective gyromagnetic factor of TM2+ ion, $`S`$ is the spin, and $`\rho `$ is the mass density calculated from the lattice parameters. The constant $`\mathrm{\Theta }_0`$ is related to the exchange integral between the nearest TM neighbors $`J_1`$,

$$\frac{2J_1}{k_B}=\frac{3\mathrm{\Theta }_0}{zS(S+1)},$$
(3)
where $`z=12`$ is the number of nearest neighbors in the wurtzite structure of Zn1-xTMxO.

A linear fit to the inverse susceptibility data intersects the $`\chi ^1=0`$ axis at a negative temperature. This result indicates the presence of antiferromagnetic interactions in the Zn1-xTMxO samples. At lower temperatures, inverse ac susceptibility deviates from the linear dependence toward a temperature close to zero. This is a result of additional antiferromagnetic interactions between the next nearest neighbor TM ions. Spalek86 . The inset to Fig. 6 shows the molar Curie constants and the Curie-Weiss temperatures for Zn1-xCoxO. From linear fits of $`C_M`$ and $`\mathrm{\Theta }`$ as a function of $`x`$, we have determined the parameters $`C_0`$ and $`\mathrm{\Theta }_0`$ for the studied TM-doped ZnO samples. The constant $`C_0`$ contains both the spin of the TM2+ ions and the effective $`g_{\mathrm{eff}}`$ factor. Mn2+-substituted ZnO was already discussed in Ref. Kolesnik02  Due to the half-filled $`3d`$ shell of Mn2+, its orbital momentum $`L=0`$, the spin is $`S=5/2`$, and the $`g_{\mathrm{eff}}`$ factor is very close to the free-spin value $`g=2`$. This is similar to other Mn-based diluted magnetic semiconductors. For Co2+, the splitting of the energy levels due to the strong spin-orbit interaction results in an increase of the effective $`g_{\mathrm{eff}}`$ factor. The values of $`g_{\mathrm{eff}}=`$2.24-2.31 were reported for Co2+ in different diluted magnetic semiconductors. Ham60 ; Koidl77 ; Lewicki89 ; Seong01 ; Alawadhi02  An enhanced $`g_{\mathrm{eff}}`$ factor is also expected for Fe2+ due to its large orbital momentum ($`L=2,S=2`$). Twardowski88

In Table I, we compare the results determined from ac susceptibility for our samples with different TM2+. Here we assume the spin $`S`$ to be equal to the theoretically expected value for each particular TM2+ ion. Our values of $`\mathrm{\Theta }_0`$ and $`2J_1/k_B`$ are similar to other TM-containing semiconductors in case of TM = Mn Spalek86 ; Furdyna88  and TM = Co. Alawadhi02  The effective $`g_{\mathrm{eff}}`$ factors for TM = Co and Fe are substantially larger than the values from literature for other diluted magnetic semiconductors doped with respective TM2+ ions. The origin of this large $`g_{\mathrm{eff}}`$ factor is not clear at the moment. The presence of a small number of TM3+ in our samples (not detectable by X-ray diffraction) is one possibility. These ions can increase the effective spin above the assumed value for TM2+. This effect is especially crucial in the case of Fe substituted ZnO which we discuss below.

The influence of synthesis conditions on the structural and magnetic properties of Zn0.95Fe0.05O is illustrated in Fig. 7. The sample annealed at 950C shows the presence of magnetite-like (Zn,Fe)3O4 impurity as displayed in the right inset. The ac susceptibility of this sample is much larger than the paramagnetic susceptibility of pure Zn0.95Fe0.05O (annealed at 1000C). The temperature dependence of ac susceptibility for the sample containing (Zn,Fe)3O4 impurity is flat at higher temperatures, which is related to a ferromagnetic order with the Curie temperature equal to 440 K. A drop of the ac susceptibility around 20 K is probably related to the low temperature Verwey transformation, similar to the one observed in stoichiometric Fe3O4 at 120 K. Bickford53  The left inset to Fig. 7 shows the magnetization curves at $`T=300`$ K for both samples. The (Zn,Fe)3O4-containing sample shows a narrow ferromagnetic hysteresis curve and the saturation magnetization close to $`0.05\mu _B`$ per formula unit, which is equal to $`1\mu _B`$ per substituting Fe ion. For ferrimagnetically coupled Fe ions in Fe3O4 one would expect the saturation magnetization $`1.33\mu _B`$ per Fe ion. This suggests that most of the Fe in the multi-phase Zn0.95Fe0.05O sample is precipitated as a (Zn,Fe)3O4 impurity. The single-phase Zn0.95Fe0.05O sample shows a linear paramagnetic dependence on the magnetic field with a very small ferromagnetic hysteretic contribution. We were able to determine the magnitude of this ferromagnetic contribution by the following numerical procedure. Linear fits of the $`M(H)`$ data were performed in the high-field ranges 40 to 70 kOe and -40 to -70 kOe, where the ferromagnetic contribution to the magnetization is expected to be fully saturated. The average of the two fitted lines was subtracted from the measured $`M(H)`$ dependence. As a result, we obtained a small ferromagnetic hysteresis curve. The saturation magnetization of this ferromagnetic contribution gives an estimate for the upper limit of the ferromagnetic (Zn,Fe)3O4 impurity of 0.1%.

## V Summary

In summary, by studying structural and magnetic properties of transition metal substituted polycrystalline ZnO samples, we have found that no bulk ferromagnetism exists for single-phase materials, contrary to the existing theories. Stoichiometric samples demonstrate paramagnetic behavior with antiferromagnetic interactions similar to other semimagnetic semiconductors substituted with transition metal ions. Selected synthesis conditions lead to impurity-induced ferromagnetic order in Co- and Fe-doped ZnO. Attempts to induce $`p`$-type doping by high-pressure oxygen synthesis lead to formation of secondary phases for which transition metals Mn, Fe, and Co exist in high oxidation states.

###### Acknowledgements.

This work was supported by NSF (DMR-0302617), the U.S. Department of Education, and the State of Illinois under HECA. The EDXS analysis was performed in the Electron Microscopy Center, Argonne National Laboratory, Argonne, IL.