text:
# Twenty Years of Timing SS433

## 1 Introduction

SS433 is the first known example of a Galactic relativistic jet source, and thus the forerunner of modern microquasar astrophysics. The optical spectrum of this object shows a number of strong, broad emission lines of the Balmer and HeI series, as well as several lines at unusual wavelengths. These latter have been identified as red/blue-shifted Balmer and HeI emission from collimated jets with intrinsic velocities of $`v0.26c`$ (Abell & Margon, 1979). Furthermore, the Doppler shifts of these features change with time in a cosinusoidal manner, leading to the label of “moving lines”. This behavior is now widely accepted to be a symptom of precession of the jet axis in SS433 on a timescale of $`164`$ days (Margon, 1984).

Early studies of the precession in SS433 indicated possible instabilities or drifts in the precessional clock (Anderson et al., 1983), which could give considerable insight into the accretion processes which must provide the precessional torque. However, Margon & Anderson (1989) reviewed ten years of SS433 timing data and concluded that while significant deviations from cosinusoidal behavior exist in SS433, the evidence for systematic long-term drifts (e.g. precessional period derivative, $`\dot{P}`$) remained inconclusive.

In this paper, we take the data set considered by Margon & Anderson (1989) and add to it more than 50 Doppler shift measurements spread over 10 years, including 9 Doppler shifts measured in 1999. Combined, these observations span more than 20 years, and thus provide an excellent data set for constraining long-term drifts in SS433’s precessional clock. We discuss the observations in Section 2. In Section 3, we present analyses of the entire data set in the context of the “kinematic model” for SS433’s precessing jets. In Section 4, we discuss the results of these analyses, and in Section 5 we present our conclusions.

## 2 Observations

The primary observations used here are optical spectroscopic observations of SS433, from a wide range of telescopes and instruments (see Margon & Anderson (1989) and references therein for details). The net result of these observations, spread over the period from June 1978 to July 1992, is the measurement of 433 Doppler shifts for the “receding jet” ($`z_1`$) and 482 Doppler shifts for the “approaching jet” ($`z_2`$) for the optical moving lines in SS433.

We obtained further spectra of SS433 in July, 1999 using the Hartung-Boothroyd Observatory (HBO) 24-inch telescope and optical CCD spectrograph. We used a 600 lines/mm grating and $`6\mathrm{}`$ slit providing a resolution of $`R800`$ ($`6\mathrm{\AA }/\mathrm{pix}`$). We present a typical spectrum in Figure 1.

We determined the Doppler shift of each HBO spectrum using only the moving $`H\alpha `$ lines, and we did so by fitting a Gaussian profile to the red and blue components separately. Note that the profiles of the moving lines are broad, time-variable, and often asymmetric. This is due to the time overlap of multiple discrete emission components, commonly referred to as “bullets” (Vermeulen et al., 1993), with typical lifetimes of $`3`$ days. These systematic deviations introduce a relatively large uncertainty in the Doppler shift determination. Based upon examination of many spectra, we find that the typical full-width at half-maximum (FWHM) is $`\mathrm{\Delta }z0.003`$, and we adopt this as our uncertainty in the Doppler shift determination $`\sigma _z`$.

## 3 Analysis

### 3.1 The Kinematic Model

Throughout our analysis of these data, we adopt the “kinematic model” for the moving lines, which assumes that the changing Doppler shifts arise from the precession of the jet axis in SS433. The simplest form of the kinematic model takes into account five components: the jet velocity $`\beta =\frac{v}{c}`$; the jet angle from the precessional axis $`\theta `$, the inclination angle of the system with respect to the observer’s line of sight $`i`$, the precession period $`P`$, and the epoch of zero precessional phase $`t_0`$. The period and zero-phase epoch combine to give the precessional phase $`\varphi =\frac{tt_0}{P}`$. The resulting Doppler shifts obey the equation

$$z_{1,2}=1\gamma [1\pm \beta \mathrm{sin}\theta \mathrm{sin}i\mathrm{cos}\varphi \pm \beta \mathrm{cos}\theta \mathrm{cos}i]$$

where $`\gamma =(1\beta ^2)^{1/2}`$. SS433 exhibits “nodding” of the jets on a $`6.5`$-day period (Katz et al., 1982) due to the $`13`$-day binary motion of SS433 (Crampton et al., 1980), which are not accounted for in this model. However, the effects of this nodding are essentially negligible for long timescale studies of the jets such as ours. We further mitigate the impact of nodding by applying a 7-day boxcar smoothing filter to the individual Doppler shifts determined above. We then used chi-squared minimization to find the best-fit parameters for the kinematic model (Table 1). The resulting model fit is plotted versus time along with the data and residuals in Figures 2-3 for $`z_1`$ and $`z_2`$. We plot the same model fit, data, and residuals versus precessional phase in Figure 4.

The resulting fit has a chi-squared residual per degree of freedom of $`\chi _\nu ^2=8.9`$, indicating the presence of statistically significant residuals. However, we can still use this fit to estimate uncertainties in the kinematic model parameters as follows. First, we scale all of the $`\sigma _z`$ values by $`\sqrt{8.9}`$, so that the residuals have $`\chi _{\nu fix}^2=1.0`$, essentially by fiat. We then take the uncertainties to be the range of a model parameter which introduces a total change of $`\mathrm{\Delta }\chi _{fix}^2=1.0`$. We also report these values in Table 1. This rescaling approach for deriving the model parameter uncertainties is statistically valid in a strict sense only if the residuals are consistent with Gaussian noise and are not correlated with any model parameters in a systamtic way. If so, then the residuals would simply indicate that we have ignored one or more sources of noise in the system when estimating the uncertainties in the individual Doppler shifts. As we show below, this is largely true, though we see some evidence of small (but statistically-significant) systematic deviations from the kinematic model. Thus, the uncertainties in the model parameters in Table 1 are likely to be good, but not perfect, statistical estimates. For the remainder of the paper, we adopt the best-fit model parameters presented in Table 1.

### 3.2 Doppler Shift Residuals

One obvious feature of Figures 2-4 is that the residuals to the model fit greatly exceed the uncertainties in the Doppler shift determinations (as also shown by the large value of $`\chi _\nu ^2`$ above). We also notice no obvious trend in the residuals versus time as would be expected for systematic timing effects, such as a constant precessional period time derivative $`\dot{P}`$. Such large, apparently random residuals have been noticed in previous timing studies of SS433 (Anderson et al. (1983); Margon & Anderson (1989)).

#### 3.2.1 Correlations in Residuals

Previous studies have also noticed that the velocity residuals in SS433 show a pattern of correlation between $`z_1`$ and $`z_2`$ (Margon & Anderson, 1989). Specifically, when we plot the residuals of $`z_1(obs)z_1(mod)`$ versus $`z_2(obs)z_2(mod)`$, we find that most of the points lie in the second and fourth quadrants (Figure 5). In other words, when the absolute value of $`z_1`$ is greater than expected, the absolute vlaue of $`z_2`$ is also greater than expected, and vice versa. The number of data points with $`z_1`$ and $`z_2`$ residuals in quadrants 2 and 4 is $`271\pm 16`$, while in quadrants 1 and 3 the number is $`110\pm 10`$ – a $`>8\sigma `$ difference.

The linear correlation coefficient between the residuals is $`r=0.69\pm 0.02`$. We estimate the uncertainty from a Monte Carlo simulation as follows. We take the 381 pairs of $`z_1`$ and $`z_2`$ residuals and add to each a random number drawn from a Gaussian distribution with mean of zero and a standard deviation of 0.003 – the typical uncertainty in the Doppler shift measurements. We then calculate the correlation coefficient of the resulting simulated distribution. We repeat this procedure 1000 times and then take the standard deviation in the correlation coefficient as the uncertainty above, $`\sigma _r=\pm 0.02`$.

This correlation pattern could have several physical sources. The effect considered most commonly in previous studies (e.g. Margon & Anderson (1989)) is that of phase noise in the precessional motion, with strict symmetry between $`z_1`$ and $`z_2`$. As the jet precessional phase either lags or leads the model ephemeris, the projected velocity amplitudes of the jets on the observers line of sight will either exceed or fall short of the model prediction. Another possible physical explanation is modulation of the velocity amplitude – “$`\beta `$-noise” – in a system which otherwise follows the 5-parameter kinematic model ideally (e.g. Milgrom et al. (1982)).

#### 3.2.2 Phase-Dependence of Residuals

Another factor which could impact the $`z_1/z_2`$ residual correlation are phase-dependent residuals to the kinematic model. Margon & Anderson (1989) found no evidence for such phase-dependence in their data. In analyzing this data set, we divided the data into 10 evenly-spaced phase intervals and calculated the average and standard deviation of the residual Doppler shifts in each bin (Table 2). None of the average residuals from the kinematic model is as large as the rms deviation for residuals in that phase bin, indicating that phase-dependent deviations from the kinematic model do not dominate the residuals. Furthermore, if we compare the average residuals of the individual phase bins to the rms scatter of all the phase bins, none of them are more than $`2\sigma `$ outliers. Thus, the amplitude of any average deviation from the model velocity does not seem to be phase-dependent.

However, when we calculate the uncertainty in the average residual for each phase bin, equal to the standard deviation in the residuals divided by the square-root of the number of points, we find that the deviations are in fact statistically significant (Table 2). That is, while the scatter around the kinematic model in any given phase bin is not dominated by systematic deviations from the model, such deviations are present in the data set. The nature of these systematic deviations are not clearly determined. Figures 2-3 show that the velocity residuals show some correlation on timescales of weeks or months. Thus, sparse sampling of the velocities combined with such correlations in the residuals could be one explanation for the apparent systematic deviations.

We used the average residuals from the kinematic model as correction factors for the quadrant analysis presented above and in Figure 5, taking the observations and subtracting both the kinematic model and the average deviation for all points within $`\mathrm{\Delta }\varphi =\pm 0.05`$ cycles of each data point. As a result, the number of $`z_1`$,$`z_2`$ residual pairs in quadrants 2 and 4 does decrease, but only to 261 (with 120 pairs in quadrants 1 and 3). Thus the correlation between $`z_1`$ and $`z_2`$ residuals remains highly statistically significant even after correction for the systematic deviations.

We also note that the relative phase-independence of the residuals raises questions regarding the nature of the residuals. As can be seen from the equation for Doppler shifts in the kinematic model, every parameter which could be time-variable ($`\beta ,P,i,\theta ,`$) either feeds directly into the phase $`\varphi `$ or is multiplied by $`\mathrm{cos}\varphi `$. Thus, noise in these terms or in $`\varphi `$ itself should result in a cosinusoidal modulation in the RMS of the Doppler shift residuals with $`\varphi `$.

#### 3.2.3 The Phase Noise Model

As mentioned above, the most commonly-invoked physical model for the velocity residuals in SS433 is “phase jitter” in the jet precession. Since the precession phase affects both jets similarly, it naturally explains the correlation between the $`z_1`$ and $`z_2`$ residuals. If such jitter can occur over timescales of weeks or months, it can also explain the long-term residual correlations evident in Figures 2-3.

We analyzed this SS433 data set following the example of Margon & Anderson (1989), determining phase errors from the velocity residuals above. We simply defined the phase error to be the phase difference between the actual phase of the observation given its epoch and the kinematic model parameters in Table 2 and the closest model point with the same observed velocity. As can be seen in Figure 4, some observed velocity amplitudes exceed the maximum model velocity amplitude, and such points were dropped from this analysis. We then divided the data set into 10-day intervals and calculated the average and standard deviation of the phase errors from all phase measurements doing that interval (including both $`z_1`$ and $`z_2`$). For 10-day intervals with only 1 phase measurement we have no estimate of the standard deviation, and thus dropped such intervals from the analysis. We plot the resulting phase noise measurements in Figure 6. We repeated this same analysis using a 30-day interval for averaging, with the results shown in Figure 7.

We note that while there are occasional trends in the residuals on timescales of several hundred days, no obvious trend is apparent over the full time span in either panel of Figure 7. The 1999 data are marginally inconsistent with zero phase residual (at the $`2.8\sigma `$ level for one of the two data points in Figure 6b). However, it is clear that this phase residual is less than many prior apprently secular deviations from the kinematic model in Figures 6-7. If we assume that some period derivative is present in SS433 over the span of our observations, these secular deviations could mask its effects up to $`\mathrm{\Delta }\varphi 0.05`$ cycles. Given the span of our observations, this corresponds to an upper limit on the period derivative of $`\dot{P}<5\times 10^5`$.

#### 3.2.4 The Velocity Amplitude Noise Model

As mentioned above, an alternate physical explanation for the velocity residuals in Figures 2-4 is noise in the intrinsic velocity of the jets. To investigate this possibility further, we calculated the intrinsic jet velocity necessary to match each observed Doppler shift, given $`\theta ,i,t_0,\mathrm{and}p`$ from the best-fit parameter set in Table 1. We plot the corresponding values for $`\beta =\frac{v}{c}`$ versus time in Figure 8 and versus precessional phase in Figure 9.

The average value of $`\beta `$ we find is 0.254 with a standard deviation of $`0.024`$. Given 507 independent measurements of $`\beta `$ in this way, we arrive at an average value of $`\beta _{ave}=0.254\pm 0.0011`$. While this value is only $`4\%`$ lower than the value given in Table 1, the difference is statistically significant at the $`7.9\sigma `$ level. This may indicate that “noise” in the Doppler shifts may in fact be impacting the parameter estimates for the kinematic model in a systematic way, as discussed in Section 3.1.

## 4 Discussion

### 4.1 Phase Noise

As noted above, the upper limit on precessional period derivative of $`\dot{P}<5\times 10^5`$ shows that there is no large long-term drift in the precessional timing properties of SS433. The presence of jitter in the system implies some “torque noise” in the process driving the precession, according to the phase noise model. However, if this were the case, that noise must average out over timescales of $`20`$ years. We can also see from Figure 7 that there are fairly large phase deviations of $`\mathrm{\Delta }\varphi 0.1`$ cycles over timescales as short as $`10`$ days. This implies that the torque noise $`\mathrm{\Delta }\tau `$ has a maximum relative amplitude of at least

$$\frac{\mathrm{\Delta }\tau _{max}}{\tau }\frac{(\mathrm{\Delta }\varphi /\mathrm{\Delta }t)}{(\varphi /t)}$$

$$1.6$$

Thus, the variation in torque can in fact exceed the time-averaged torque driving the precession. This may be a problem for certain physical models of the precession and timing noise in SS433.

Finally, we note that the phase noise model is incapable of producing the observed Doppler shifts which exceed the maximum amplitude predicted by the kinematic model. We have considered the possibility that the phase noise itself causes the $`\chi `$-squared fitting procedure used to determine the model parameters to systematically underestimate the true jet velocity, and thus “undershoot” the maxima. However, Monte Carlo simulations of data sets with higher true velocities and phase noise identical to that observed here fail to produce such undershooting. Therefore, we conclude that phase noise model cannot reproduce the observed Doppler shift residuals near the maximum projected velocities.

### 4.2 Velocity Noise

The alternate “$`\beta `$-noise” model, on the other hand, can clearly explain the excess velocity at the extrema (and any other precessional phase) by the changing jet velocity amplitude. Such a model also has a physical basis, given recent advances in the modeling of relativistic jet production. Meier et al. (2001) dicuss a scenario where such jets are launched by a magnetic accretion disk instability around a black hole (or other compact object). Variations in the accretion flow onto the compact object (i.e. $`\dot{M}`$, intrinsic magnetic field, etc.) can alter the radius at which the magnetic field saturates and the jet is launched, and thus the jet velocity. The relation between jet velocity and launch radius for a non-rotating black hole follows:

$$\beta (R)=\sqrt{\frac{2R_g}{R}}$$

where $`\beta =\frac{v}{c}`$, and $`R_g`$ is the gravitational radius of the black hole (one-half of the Schwarzschild radius).

### 4.3 Jitter models and phase-dependence of residuals

It is also interesting to view these model in light of the apparent lack of phase-dependence in the Doppler shift residuals noted above. By differentiating the equation for Doppler shifts in the kinematic model (eqn. 1) with respect to the potentially time-varying model components ($`\beta `$, $`i`$, $`\theta `$, and $`\varphi `$) we can see the relative phase-dependence of Doppler shift residuals on deviations in each term. In the phase noise model, we would expect the following dependence:

$$\mathrm{\Delta }z=(\gamma \beta \mathrm{sin}\theta \mathrm{sin}i\mathrm{sin}\varphi )\mathrm{\Delta }\varphi $$

$$0.1\mathrm{sin}\varphi \mathrm{\Delta }\varphi $$

Thus, we would expect the amplitude of the Doppler shift residuals to be sinusoidally modulated with phase, in apparent contradiction with our analyses above.

For the “$`\beta `$-noise” model, we have:

$$\mathrm{\Delta }z=(\gamma \mathrm{sin}\theta \mathrm{sin}i\mathrm{cos}\varphi \gamma \mathrm{cos}\theta \mathrm{cos}i)\mathrm{\Delta }\beta $$

$$(0.35\mathrm{cos}\varphi 0.2)\mathrm{\Delta }\beta $$

(with the approximation that $`\gamma /\beta 0`$). Again, we have a modulation of the Doppler shift residual amplitude dominated by a term varying cosinusoidally with respect to phase.

For variations in the angle between the jet axis and the precessional axis, $`\theta `$, we have:

$$\mathrm{\Delta }z=(\gamma \beta \mathrm{cos}\theta \mathrm{sin}i\mathrm{cos}\varphi +\gamma \beta \mathrm{sin}\theta \mathrm{cos}i)\mathrm{\Delta }\theta $$

$$(0.24\mathrm{cos}\varphi +0.09)\mathrm{\Delta }\theta $$

again dominated by a $`\mathrm{cos}\varphi `$ term.

Finally, for variations in the system inclination angle, $`i`$, we have:

$$\mathrm{\Delta }z=(\gamma \beta \mathrm{sin}\theta \mathrm{cos}i\mathrm{cos}\varphi +\gamma \beta \mathrm{cos}\theta \mathrm{sin}i)\mathrm{\Delta }i$$

$$(0.02\mathrm{cos}\varphi +0.24)\mathrm{\Delta }i$$

Interestingly, the Doppler shift residual amplitudes for “$`i`$-noise” would be dominated by a constant term, with only a small dependence on phase. Thus, this is the only parameter in the kinematic model for which variations causing the Doppler shift residuals are consistent with their observed phase-independence. Unfortunately, we know of no physical model for such variability at this time. Furthermore, based on the mass estimates for the compact object and companion star ($`10M_{}`$ total), the known 13.5-day binary period, and the presence of deviations as large as 0.01 rad/day in inclination angle, the change in rotational energy would require average powers of $`\dot{E}10^3L_{Edd}`$, which seems implausible.

One other possible explanation is that the Doppler shift residuals are due to variations in both precessional phase and one (or more) of the other model parameters. In that case, the residual amplitude would depend on a sum of $`\mathrm{cos}\varphi `$ and $`\mathrm{sin}\varphi `$ terms which could potentially smooth out any phase-dependence in the residuals.

## 5 Conclusions

We have presented observations of the Doppler-shifted optical moving lines in SS433 spanning over 20 years. We draw the following conclusions based on the data:

* We find parameters for the kinematic model for the jet precession which are similar to those found by previous authors (e.g. Margon & Anderson (1989)).
* We find a strong correlation between residuals to the models fits for the two jets $`z_1`$ and $`z_2`$, with a linear correlation coefficient of $`r=0.69\pm 0.02`$.
* We find that the residuals to the kinematic model fit are not dominated by systematic phase-dependent deviations from the model. However, systematic phase-dependent deviations from the kinematic model are seen in the data set at a low level.
* If we adopt a “phase noise” model for the velocity residuals, we find correlated deviations over timescales of months to years, but no long-term trend over the full data set. We place a limit on the precessional period derivative of $`\dot{P}<5\times 10^5`$.
* Noise in any single parameter of the kinematic model seems unable to explain the observed phase-independence of the velocity residuals in SS433. However, variations in both phase and one of the other parameters would vary as the weighted sum of $`\mathrm{cos}\varphi `$ and $`\mathrm{sin}\varphi `$ terms, which could smooth out any phase-dependence of the Doppler shift residuals.

The authors thank D. Lai and I. Wasserman for helpful discussions of these results. SSE is supported in part at Cornell by an NSF CAREER award (NSF-9983830), and PBC was partially supported by this grant. DMK was supported at Cornell by a NASA Space Grant summer research fellowship.


format: text

text:
# A note on cluster methods for strongly correlated electron systems

\[

## Abstract

We develop, clarify and test various aspects of cluster methods dynamical mean field methods using a soluble toy model as a benchmark. We find that the Cellular Dynamical Mean Field Theory (C-DMFT) converges very rapidly and compare its convergence properties with those of the Dynamical Cluster Approximation (DCA). We propose and test improved estimators for the lattice self energy within C-DMFT.

\]

The development of dynamical mean field methods has resulted in significant advances in our understanding of strongly correlated electron systems, in particular in the area of the Mott transition . This method, captures the local effects of correlations such as the Kondo effect and the transfer of spectral weight between the coherent and the incoherent part of the spectral function. It suffers however from limitations arising from its single site mean field character such as the lack of k dependence of the self energy. In the context of disordered system, the dynamical mean field theory has a precursor in the famous coherent potential approximation (CPA) for which cluster extensions such as the molecular CPA  have been formulated. Naturally, generalizations of these and other statistical mechanical approaches such as the Bethe Peirls approximation to the area of quantum interacting systems, have been investigated recently , ,  , ,, . This area of investigation is in its beginning stages, and comparative studies of the various methods are important to increase our understanding of the strengths and the limitations of these methods, at a level comparable to our present understanding of the single site dynamical mean field theory. This note is a contribution in this direction. We focus, here on the CDMFT  and the DCA  method, because both methods have been proved to be manifestly causal, i.e. the output of an approximated solution of the cluster equations is causal, as long as a causal method is used for the solution of the impurity model. Because no numerical investigations of the C-DMFT have been carried out, we test the performance of this method in a simple soluble model that was introduced by I. Affleck and B. Marston . It has a k dependent, albeit static, self energy, and therefore is a simple playground to explore the cluster method.

In the first part of the paper we describe the CDMFT , and introduce a real space formulation of the DCA. to facilitate the comparison with CDMFT. In the second part we compare the predictions of DCA and CDMFT for the short distance behavior of correlation functions, in different cluster sizes against the exact solution. In the CDMFT approach the lattice self energy is a derived quantity which needs to be estimated from the cluster self energy, an auxiliary quantity which enters the dynamical mean field equations. In the second part of the paper we provide improved estimators for the lattice self energy and discuss how they improve the convergence to the exact answer as a function of the cluster size.

Real space formulation the cluster schemes: A fairly general model of strongly correlated electrons contains hopping and interaction terms. It is defined on a lattice of $`L^d`$ sites in d dimensions, and we divide the lattice in $`(L/L_c)^d`$ cubic clusters of $`L_c^d`$ sites (more general forms can also be considered). We denote with $`𝐞_i`$ the internal cluster position and with $`𝐑_n`$ the cluster position in the lattice (therefore the position of the i-th site of the n-th cluster is $`𝐑_n+𝐞_i`$). The lattice Hamiltonian is expressed in terms of fermionic operators $`f_{𝐑_n,\alpha }^{}`$ and $`f_{𝐑_m,\beta }`$ and can be written as:

$`H`$ $`=`$ $`{\displaystyle \underset{n,\alpha ,m,\beta }{}}t_{\alpha ,\beta }(𝐑_n𝐑_m)f_{𝐑_n,\alpha }^{}f_{𝐑_m,\beta }`$ (1)
$`+`$ $`{\displaystyle \underset{n,\alpha ,m,\beta ,n^{},\alpha ^{},m^{},\beta ^{}}{}}U_{\alpha ,\beta ,\alpha ^{},\beta ^{}}(\{𝐑\})f_{𝐑_n,\alpha }^{}f_{𝐑_m,\beta }f_{𝐑_n^{},\alpha ^{}}^{}f_{𝐑_m^{},\beta ^{}}`$ (2)
$`\alpha =i,\sigma `$ and $`\sigma `$ is an internal degree of freedom (i.e. a spin, spin-orbital or band index). Most cluster schemes, neglect the interaction terms between different clusters. The effects of those interactions, can be treated using the extended dynamical mean field approach  but we will not discuss these improvements in this paper. All the different cluster schemes can be formulated as a self consistent equation for the cluster self-energy which consist of the following loop:(i) start with a guess of the cluster self-energy $`(\mathrm{\Sigma }_c)_{\alpha ,\beta }`$, (ii) from the cluster self energy compute the Weiss function or host cluster propagator $`(G_0)_{\alpha ,\beta }`$, which enters in the effective action for the cluster degrees of freedom, (iii) use the effective action compute the cluster Green function $`(G_c)_{\alpha ,\beta }`$, (iiii) compute the new cluster self-energy, (iiiii) iterate this loop until the convergence is reached. The DCA and CDMFT schemes differ in the way step (ii) is carried out. Within the CDMFT one obtains the Weiss function from the cluster self energy by the equation:

$$\widehat{G}_0^1=\left(\left(\frac{L_c}{L}\right)^d\underset{𝐊}{}\frac{1}{(i\omega _n+\mu )\widehat{}\widehat{t}(𝐊)\widehat{\mathrm{\Sigma }}_c}\right)^1+\widehat{\mathrm{\Sigma }}_c$$
(3)
where $`t(𝐊)_{\alpha ,\beta }`$ is the Fourier transform of the hopping matrix in eq. (1) with respect to $`𝐑_n𝐑_m`$ ($`𝐤`$ is a wave-vector in the Brillouin zone reduced by $`L_c`$ in each direction), $`\omega _n`$ is the Matsubara frequency and $`\mu `$ is the chemical potential. Once the Weiss function has been computed one can obtain $`\widehat{G}_c`$ by functional integration of the single site action. Step (iiii) is carried out using the definition of the cluster self-energy: $`\widehat{\mathrm{\Sigma }}_c=\widehat{G}_0^1\widehat{G}_c^1`$.

To facilitate the comparison with the CDMFT, in the following we shall derive the DCA scheme using the real space formulation of the cluster. To lighten the notation we will assume that the variable $`\sigma `$ is conserved to make all the cluster matrices diagonal in $`\sigma `$ and subsequently we will drop this index. We take periodic boundary condition on the cluster and we define the matrix $`U_{i,j}(𝐊)=\mathrm{exp}(i𝐊𝐞_i)\delta _{𝐞_i,𝐞_j}`$. It is crucial for the following to note that the matrix $`\widehat{t}(𝐊)`$ has the following representation: $`t(𝐊)_{i,j}=\frac{1}{L_c^d}_{𝐤_c}e^{i(𝐊+𝐤_c)(𝐞_i𝐞_j)}t(𝐊+𝐤_c)`$, where $`𝐤_c`$ are the cluster momenta. Therefore the matrix $`\widehat{U}(𝐊)\widehat{t}(𝐊)\widehat{U}^{}(𝐊)`$ is diagonal with respect to cluster momenta: $`(\widehat{U}(𝐊)\widehat{t}(𝐊)\widehat{U}^{}(𝐊))_{i,j}=\frac{1}{L_c^d}_{𝐤_c}e^{i𝐤_c(𝐞_i𝐞_j)}t(𝐊+𝐤_c)`$. Using this property one can write the (ii) DCA eq. in real space as:

$$\widehat{G}_0^1=\left(\left(\frac{L_c}{L}\right)^d\underset{𝐊}{}\frac{1}{(i\omega _n+\mu )\widehat{}\widehat{U}\widehat{t}\widehat{U}^{}(𝐊)\widehat{\mathrm{\Sigma }}_c}\right)^1+\widehat{\mathrm{\Sigma }}_c$$
(4)
Since the matrices $`\widehat{U}\widehat{t}\widehat{U}^{}(K)`$ and $`\widehat{\mathrm{\Sigma }}_c`$ are diagonal with respect to $`𝐤_c`$, this equation coincides with the DCA equations of Jarell et. al.  after a Fourier transformation with cluster momenta. Once $`\widehat{G}_0`$ is known, $`\widehat{G}_c`$ is computed by functional integration of the cluster effective action and the new cluster self-energy is obtained by $`\mathrm{\Sigma }_c(𝐤_c)=G_0^1(𝐤_c)G_c^1(𝐤_c)`$.
Eq. 4 allows a direct formulation of DCA in real space and a detailed comparison with CDMFT. We also note this real space formulation can be used to defined many causal cluster schemes, by introducing a different matrix $`U(k)`$ in the previous equation.

A simplified one dimensional large-N model: comparison between the exact solution and the predictions of the cluster schemes. In the following we focus on a simple one dimensional model, originally introduced and studied by I. Affleck and B. Marston  in two dimensions. We compare the DCA and CDMFT schemes to its exact solution. This model is a generalization of the Hubbard-Heisenberg model where the SU(2) spins are replaced by a SU(N) spins, the on site repulsion is scaled as $`1/N`$ and the large N limit is taken. Its Hamiltonian reads:

$`H`$ $`=`$ $`t{\displaystyle \underset{i,\sigma }{}}(f_{i,\sigma }^{}f_{i+1,\sigma }+f_{i+1,\sigma }^{}f_{i,\sigma })`$ (5)
$`+`$ $`{\displaystyle \frac{J}{2N}}{\displaystyle \underset{i,\sigma ,\sigma ^{}}{}}(f_{i,\sigma }^{}f_{i,\sigma ^{}}f_{i+1,\sigma ^{}}^{}f_{i+1,\sigma }+f_{i+1,\sigma }^{}f_{i+1,\sigma ^{}}f_{i,\sigma ^{}}^{}f_{i,\sigma })`$ (6)
where $`i=1,\mathrm{},L`$ and $`\sigma =1,\mathrm{},N`$ and we take the large $`L`$ and $`N`$ limits. In the following we will use $`J`$ as the unite of temperature and therefore we put $`J=1`$ and we rescale the hopping term $`tt/J`$. The thermodynamics of this model can be solved exactly since in the large $`N`$ limit the quantity $`\chi =\frac{1}{N}_\sigma f_{i,\sigma }^{}(t)f_{i+1,\sigma }(t)`$ does not fluctuate. Indeed (5) reduces to a free-fermions Hamiltonian with a “renormalized” hopping term $`tt+\chi `$ and a self-consistent condition on $`\chi `$:

$$\chi =\frac{1}{L}\underset{k}{}f(\beta E_k)\mathrm{cos}k,E_k=2(t+\chi )\mathrm{cos}k+\mu $$
(7)
where $`\mu `$ is the chemical potential, $`f(\beta E_k)`$ is the Fermi function and $`\beta `$ is the inverse temperature.
We now apply the DCA approximation to the Hamiltonian (5). As previously, the computation is simplified by the fact that the quantity $`\chi _{cl}^{DCA}=\frac{1}{N}_\sigma f_{i,\sigma }^{}(t)f_{i+1,\sigma }(t)`$, where $`i`$ and $`i+1`$ belong to the same cluster, does not fluctuate in the large-N limit. As a consequence the functional integral on the cluster degrees of freedom reduce to a simple Gaussian integral. Thus, imposing periodic boundary condition on the cluster, the eq. $`𝒢^1(k_c)=G_c^1(k_c)+\mathrm{\Sigma }_c(k_c)`$ implies: $`\mathrm{\Sigma }_c(k_c)=2\chi _{cl}^{DCA}\mathrm{cos}k_c`$. Using the second DCA equation which expresses the cluster Green-function as a function of $`\mathrm{\Sigma }_c(k_c)`$, we obtain the self-consistent DCA relation for $`\chi _{cl}^{DCA}`$:

$`\chi _{cl}^{DCA}`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{K,k_c}{}}f(\beta E_{K,k_c})\mathrm{cos}k_c,`$ (9)
$`E_{K,k_c}=2t\mathrm{cos}(k_c+K)2\chi _{cl}^{DCA}\mathrm{cos}k_c+\mu `$
Note that in the infinite cluster limit one recovers the exact equation (7).

Now we focus on the CDMFT approximate solution. As in the DCA case, since the quantity $`(\chi _{cl}^{CDMFT})_i=\frac{1}{N}_\sigma f_{i,\sigma }^{}(t)f_{i+1,\sigma }(t)`$ does not fluctuate, one obtains: $`(\mathrm{\Sigma }_c)_{i,j}=(\chi _{cl}^{CDMFT})_i\delta _{i,j1}+(\chi _{cl}^{CDMFT})_j\delta _{i,j+1}`$. This is the generalization of the corresponding DCA expression to a case without periodic boundary condition. Note that now the quantity $`(\chi _{cl}^{CDMFT})_i`$ may depend on the cluster index. Denoting the eigenvectors and the eigenvalues of the matrix $`\widehat{t}+\widehat{\mathrm{\Sigma }}`$ respectively $`\psi _i^\nu (K)`$ and $`\lambda ^\nu (K)`$ ($`\nu =1,\mathrm{},L_c`$), the (ii) CDMFT equation, which expresses the cluster Green-function in terms of the cluster self-energy, reads:

$$(G_c)_{i,j}=\frac{L_c}{L}\underset{K,\nu }{}\psi _i^\nu (K)(\psi _j^\nu (K))^{}\frac{1}{i\omega _n+\mu \lambda ^\nu (K)}$$
(10)
Using this expression we finally get the self-consistent CDMFT equation on $`(\chi _{cl}^{CDMFT})_i`$:

$`(\chi _{cl}^{CDMFT})_i`$ $`=`$ $`{\displaystyle \frac{L_c}{L}}{\displaystyle \underset{K,\nu }{}}\psi _i^\nu (K)(\psi _{i+1}^\nu (K))^{}f(\beta E_{K,\nu }),`$ (12)
$`E_{K,\nu }=\mu \lambda ^\nu (K)`$
Notice that eq. 12 corresponds to the exact solution of a model defined by a Hamiltonian similar to (5) in which $`J_{i,i+1}`$ equals $`1`$ if $`i`$ and $`i+1`$ belong to the same cluster and zero otherwise. This implies in particular that in the infinite cluster limit the CDMFT approximation gives back the exact solution. We have numerically solved the self consistent equations (7,9,12) to compare the DCA and CDMFT predictions for different cluster sizes to the exact solution.

In Fig. (1) we plot the result of this analysis for $`t=1,\mu =1`$ as a function of $`\beta `$ and for different cluster sizes. The two methods converge (the convergence is not uniform in $`\beta `$) toward the exact solution for high enough $`L_c`$ but C-DMFT converges better that DCA. Indeed the C-DMFT result are already surprisingly good for $`L_c=2`$. However, as we shall discuss below, there are two different ways to compute $`\chi `$ within the cluster methods. The one used here is based on a real space cluster intuition. The second one, relies on a momentum space intuition and computes the correlation functions from the k dependent lattice Greens function, this procedure was advocated in ref  and . we will show that it is also gives accurate results for small clusters.

The lattice self-energy. We now address the computation of the lattice self energy. In DCA a discretized form of the lattice self energy in momentum space, enters directly in the evaluation of $`G_0`$. On the other hand, CDMFT focuses on estimating the cluster greens function, and the lattice self energy does not participate in the mean field equations, and has to be estimated later from the cluster self energy. For the simplified large-N one dimensional model studied in this paper the DCA prediction for the lattice self-energy reads:

$$\mathrm{\Sigma }_{latt}^{DCA}(k)=\mathrm{\Sigma }(k_c)=2\chi _{cl}^{DCA}\mathrm{cos}k_c$$
(13)
where $`k`$ belongs to $`[\frac{\pi }{L_c}+k_c,k_c+\frac{\pi }{L_c}]`$. Whereas for C-DMFT an estimator for the lattice self energy is constructed using the matrix $`S_{R_n,\alpha ;i}=\delta _{R_n+\alpha ,i}`$ where $`\alpha `$ is the internal cluster index and $`i`$ denotes a lattice site. The simplest form  is: $`\mathrm{\Sigma }_{latt}^{CDMFT}(k)=_{\alpha ,\beta }\stackrel{~}{S}_\alpha ^{}(k)\mathrm{\Sigma }_{\alpha ,\beta }\stackrel{~}{S}_\beta (k)`$, where $`\stackrel{~}{S}_\alpha (k)`$ is the Fourier transform of the matrix $`S`$ with respect to the original lattice index $`i`$. $`\stackrel{~}{S}_\alpha (k)`$ can be easily written in terms of the matrix $`\widehat{U}(k)`$ defined before: $`\stackrel{~}{S}(k)_\alpha =U_{\alpha ,\alpha }(k)/\sqrt{L_c}`$. Therefore the relationship between the lattice and the cluster self energy reads: $`\mathrm{\Sigma }_{latt}^{CDMFT}(k)=_{\alpha ,\beta }[U^{}(k)\mathrm{\Sigma }_cU(k)]_{\alpha ,\beta }/L_c`$. For example, in the case of the two site cluster we find: $`\mathrm{\Sigma }_{latt}^{CDMFT}(k)=\chi _{cl}^{CDMFT}cos(k)`$, whereas the exact solution gives $`\mathrm{\Sigma }_{latt}(k)=2\chi _{ex}cos(k)`$. As a consequence, even if the value of $`\chi `$ is well predicted by the CDMFT there is a factor 2 between the two self-energies. The reason of this discrepancy may be understood writing the simple estimator of the lattice self energy  in real space $`(\mathrm{\Sigma }_{latt})_{ij}=_{\alpha ,\beta :\alpha \beta =ij}(\mathrm{\Sigma }_c)_{\alpha ,\beta }/L_c`$. This means that the lattice self energy for a certain value of $`ij`$ is obtained averaging over all the cluster self energy elements corresponding to $`\alpha \beta =ij`$. In the limit of an infinite cluster translation invariance implies that the cluster self energy coincides with the lattice self energy in the bulk. Therefore the factor $`1/L_c`$ cancels and we get the exact solution. However for a finite lattice there are only $`L_c1`$ factors for $`ij=1`$, $`L_c2`$ factors for $`ij`$=2 , …$`L_ck`$ factors for $`ij=k`$. Therefore it is highly desirable to have improved estimators for smaller size clusters in which the formula in which the average over all the factors having $`\alpha \beta =k`$ is weighted by their number $`1/(L_ck)`$. One could also think to put an extra weight to extract the lattice self-energy only from the sites in the bulk, for which the CDMFT result should be better. We propose new general class of estimators for the lattice self energy in terms of the cluster self energy, that inherit its causality property.

$$(\mathrm{\Sigma }_{latt})_{ij}=\underset{\alpha ,\beta :\alpha \beta =ij}{}w_{\alpha ,\beta }(\mathrm{\Sigma }_c)_{\alpha ,\beta }$$
(14)
where the matrix $`w_{\alpha ,\beta }`$ is positive definite and $`_{\alpha ,\beta :\alpha \beta =ij}w_{\alpha ,\beta }1`$ for $`L_c\mathrm{}`$ (this guarantees a good behavior in the infinite cluster limit). Using that the trace of the product of two positive definite matrices is positive, one can easily prove that if the cluster self-energy is causal this formula produce a lattice self-energy which is also causal. Note that (14) does not change the behavior for an infinite cluster, but can really improve the results for finite cluster sizes. For example, in Fig. 2 we compare the exact lattice self-energy to the DCA and the CDMFT predictions (using the initial estimator proposed in  and the simple improvement $`w_{\alpha ,\beta }=1/(L_c1)`$ which weights in the right way at least the terms with $`\alpha \beta =1`$) for $`\beta =\mu =t=1`$. We remark that there is an excellent agreement between CDMFT and the exact solution after that our simple improvement has been taken into account already for $`L_c=2`$. Whereas the DCA prediction becomes good for $`L_c5`$.

Relation between lattice and cluster observables. Once the lattice self-energy has been obtained within a cluster method, the lattice Green function can be straightforwardly computed. This offers a different way of estimating the first neighbor correlation function $`\chi `$, using the lattice Greens function. This quantity can be computed inside the cluster ($`\chi _{cl}`$) or using the lattice Green function, obtained by the lattice self-energy, ($`\chi _{la}`$) and the two results do not coincide in general. In the case of C-DMFT one can understand what are the approximations responsible for this difference and why they are small. The C-DMFT approach is based on the cavity procedure  which, if it was carried out exactly, it would give back the same answer for the lattice and cluster observables. However, in the approximated cavity procedure adopted in the C-DMFT, one assumes that the contribution to the effective action coming from tracing out all the degrees of freedom outside the cluster is purely gaussian. This is clearly not the case in general and it is the main reason for the non-coincidence of lattice and cluster observables.

In Fig. (3) we compare the DCA and CDMFT predictions for the lattice and cluster values of $`\chi `$ to the exact solution for a two site cluster, for $`t=1,\mu =1`$ as a function of $`\beta `$. These curves display the typical behavior found also for other values of the control parameters: $`\chi _{la}^{DCA}`$ is quite better than its cluster counterpart, whereas the C-DMFT prediction is quite stable. This is probably the result of having an approximate cavity construction for the C-DMFT . Moreover we remark that $`\chi _{la}^{CMDT2}`$, obtained using the first improvement for the self-energy discussed above, almost coincides with the exact solution. Comparing the C-DMFT and the DCA lattice values of $`\chi `$ we note that C-DMFT gives usually a little bit better answer in terms of accuracy and convergence with respect to cluster size. This is perhaps due to the smoothness of the self-energy in the C-DMFT case.

It would be nice to eliminate the cluster self energy altogether from the C-DMFT approach or to use a self energy without discontinuities in DCA, in the spirit of the work of Katsnelson and Lichtenstein . However, we were unable to prove manifest causality of this approach.

In summary, in this short note we compared the performance of the DCA method with that of the cellular DMFT, in a very simple toy model. We have also proposed new estimators for the lattice self energy within CDMFT, which is more efficient. Our study shows that a direct application of CDMFT, i.e. without exploiting the flexibility inherent in the choice of basis in its most general formulation is very efficient in converging to the correct solution already for a two site cluster. Comparing the DCA and the C-DMFT predictions are a little bit better in terms of accuracy and convergence with respect to cluster size. DCA estimates of physical quantities, are most accurately carried out using the lattice Greens function, and not from the real space cluster correlation functions. This is stressed in ref , who view DCA as a momentum space method. C-DMFT is not so sensitive to the choice of cluster or lattice estimators, because of the underlying a cavity construction of the derivation . These results are very encouraging, and warrant further applications of these methods to more realistic and difficult problems. Since the most glaring deficiency of the CDMFT method is that it does not attempt to take into account in a direct fashion the translation invariance of the problem, we concentrated on the phase of this model (finite temperatures) which is translationally invariant. We expect that CDMFT will perform even better for the ground state properties since in this case translation invariance is broken by dimerization.

Acknowledgments: This work was supported by the NSF under grant DMR-0096462 and by the Rutgers Center for Materials Theory. GK acknowledges useful discussions with A. Georges, M. Jarrell, H. Krishnamurthy and Sasha Lichtenstein.


text:
# Identified Particle Elliptic Flow in Au+Au Collisions at √𝑠__{𝑁𝑁}=130 GeV

## Abstract

We report first results on elliptic flow of identified particles at mid-rapidity in Au+Au collisions at $`\sqrt{s_{_{NN}}}=130`$ GeV using the STAR TPC at RHIC. The elliptic flow as a function of transverse momentum and centrality differs significantly for particles of different masses. This dependence can be accounted for in hydrodynamic models, indicating that the system created shows a behavior consistent with collective hydrodynamical flow. The fit to the data with a simple model gives information on the temperature and flow velocities at freeze-out.

The goal of the ultra-relativistic nuclear collision program is the creation of a system of deconfined quarks and gluons qm99 . The azimuthal anisotropy of the transverse momentum distribution for non-central collisions is thought to be sensitive to the early evolution of the system sorge . The second Fourier coefficient of this anisotropy, $`v_2`$, is called elliptic flow. It is an important observable since it is sensitive to the rescattering of the constituents in the created hot and dense matter. This rescattering converts the initial spatial anisotropy, due to the almond shape of the overlap region of non-central collisions, into momentum anisotropy. The spatial anisotropy is largest early in the evolution of the collision, but as the system expands and becomes more spherical, this driving force quenches itself. Therefore, the magnitude of the observed elliptic flow should reflect the extent of the rescattering at relatively early time sorge .

Elliptic flow in ultra-relativistic nuclear collisions was first discussed in Ref. olli92  and has been studied intensively in recent years at the Alternating Gradient Synchrotron (AGS) e877flow2 ; e895 , Super Proton Synchrotron (SPS) na49prl ; na49flow ; wa98  and at the Relativistic Heavy Ion Collider (RHIC) starflow  energies. The studies at top AGS energy and SPS have found that elliptic flow at these energies is in the plane defined by the beam direction and the impact parameter, $`v_2>0`$, as expected from most models. The pion elliptic flow for relatively peripheral collisions increases with beam energy olli98  from about 0.02 at the top AGS energy e877flow2 , 0.035 at the SPS na49flow  to about 0.06 at RHIC energies starflow . The increased magnitude of the integrated elliptic flow at RHIC reaches the values predicted by hydrodynamical models, which are based on the assumption of complete local thermalization.

The first elliptic flow results from RHIC were for charged particles. The differential charged particle flow, $`v_2`$($`p_t`$), shows an almost linear rise with transverse momentum, $`p_t`$, up to 1.5 GeV/$`c`$. At $`p_t>`$ 1.5 GeV/$`c`$, the $`v_2`$($`p_t`$) values start to saturate, which might indicate the onset of hard processes starflow ; wang ; gyulassy ; molnar . The behavior of $`v_2`$($`p_t`$) up to 1.5 GeV/$`c`$ is consistent with a hydrodynamic picture. Hydrodynamics assumes complete local thermalization at the formation of the system, followed by an evolution governed by an Equation Of State (EOS). However, the $`p_t`$-integrated elliptic flow, $`v_2`$, as a function of centrality, and the differential $`v_2`$($`p_t`$), show little sensitivity to the EOS used kolb1 . Studies of the mass dependences of elliptic flow for particles with $`p_t<`$ 1.5 GeV/$`c`$ provide important additional tests of the hydrodynamical model pasi2 . Similar to the identified single particle spectra, where the transverse flow velocity has been extracted from the mass dependence of the slope parameter nxu , the $`v_2`$($`p_t`$) for different mass particles allows the extraction of the elliptic component of the flow velocity eos ; voloshin\_e877 . Moreover, the details of the dependence of elliptic flow on particle mass and transverse momentum are sensitive to the temperature, transverse flow velocity, its azimuthal variation, and source deformation at freeze-out. In this Letter we report the first results for identified particle $`v_2`$($`m,p_t`$) in Au+Au collisions at RHIC at $`\sqrt{s_{_{NN}}}=130`$ GeV.

The Solenoidal Tracker At RHIC (STAR) STAR , is ideally suited for measuring elliptic flow due to its azimuthal symmetry and large coverage. The detector consists of several sub-systems in a large solenoidal magnet. The Time Projection Chamber (TPC) covers the pseudorapidity range $`|\eta |<1.8`$ for collisions in the center of the TPC. The magnet was operated at a 0.25 Tesla field, allowing tracking of particles with $`p_t>75`$ MeV/c. Two Zero Degree Calorimeters ZDC  located at $`|\mathrm{sin}\theta |<`$ 0.002, which mainly detect fragmentation neutrons, are used in coincidence for the trigger. The TPC is surrounded by a scintillator barrel which measures the charged particle multiplicity within $`|\eta |<1`$, for triggering purposes.

For this analysis, 120,000 minimum-bias events were selected with a primary vertex position within 75 cm longitudinally of the TPC center and within 1 cm radially of the beam line. For determination of the event plane starflow , charged particle tracks were selected with $`0.1<p_t2.0`$ GeV/c. All tracks used in this analysis passed within 2 cm of the primary vertex and had at least 15 measured space points. Also, the ratio of the number of measured space points to the expected maximum number of space points for that particular track was required to be greater than 0.52, suppressing split tracks from being counted more than once. The tracks used for the determination of the reaction plane were within $`|\eta |<1.0`$, and the tracks used to calculate the elliptic flow were within $`|\eta |<1.3`$. These cuts are similar to the ones used in Ref. starflow , and the analysis results presented here are not sensitive to those cuts.

The pions, protons, and antiprotons were selected according to their specific energy loss ($`dE/dx`$) in the TPC in the transverse momentum range of 0.175 – 0.75 GeV/$`c`$, 0.5 – 0.9 GeV/$`c`$ and 0.3 – 0.9 GeV/$`c`$, respectively. The lower $`p_t`$ values were chosen such that the energy loss in the detector was negligible. At low momentum, the proton background due to secondary interactions with the detector material is significant. Therefore, only protons above a transverse momentum of 0.5 GeV/$`c`$ were used in this analysis. At a momentum ($`p`$) of 0.5 GeV/$`c`$, the $`dE/dx`$ resolution was found to be of the order of 11% for a typical long track in the STAR TPC. The kaons, because of their overlap in $`dE/dx`$ with the electrons and positrons, were selected in the momentum ranges of 0.30 – 0.40 GeV/$`c`$ and 0.60 – 0.70 GeV/$`c`$. The raw yields of the pions, kaons, protons, and antiprotons were obtained from fitting the $`dE/dx`$ distributions for each $`\eta `$, $`p`$ and centrality bin with multiple Gaussians, and requiring greater than 90% purity (particle’s yield divided by the sum of yields for all particles at that same $`dE/dx`$). Figure 1 shows the $`dE/dx`$ distributions for negative tracks in the TPC for central events and $`|\eta |0.1`$. The upper panel shows the $`dE/dx`$ distribution for the momentum range of 0.38 – 0.40 GeV/$`c`$ and the multiple Gaussian fit. In this momentum range, the $`\pi ^{}`$, $`e^{}`$, K- and $`\overline{p}`$ are clearly separated. The lower panel shows the $`dE/dx`$ distribution for the momentum range of 0.70 – 0.72 GeV/$`c`$. For the higher momentum bin, the $`\pi ^{}`$, $`e^{}`$ and K- overlap; however, the $`\overline{p}`$ is still clearly separated. Even though it is still possible to extract kaons with 90% purity in this momentum bin, we only used the kaons up to 0.70 GeV/$`c`$.

The flow analysis method meth  involves the calculation of the event plane angle, which is an experimental estimator of the real reaction-plane angle. For the analysis presented in this Letter, each particle was correlated with the event plane from all the other particles (for the other methods, see QMsnellings ).

The differential elliptic flow, $`v_2`$, depends on mass, rapidity ($`y`$) and $`p_t`$. In Fig. 2, $`v_2`$($`p_t`$) is shown for pions, kaons, and protons + antiprotons for minimum-bias collisions footminbias , integrated over rapidity and centrality by taking the multiplicity-weighted average. The uncertainties shown are statistical only. Using the same procedure to estimate the systematic uncertainties as in Ref. starflow  for $`v_2`$ integrated over $`p_t`$, the systematic uncertainty for minimum-bias data is 13%. We have verified that the positive and negative identified particles used in this analysis have the same $`v_2`$($`p_t`$) within statistical uncertainties. For $`v_2`$($`p_t`$), the pions were integrated over $`|y|1.0`$, the kaons over $`|y|0.8`$, and the protons + antiprotons over $`|y|0.5`$ (the rapidity ranges chosen correspond to approximately consistent $`|\eta |`$ coverage for all of the selected particles). Mathematically the $`v_2`$ value at $`p_t=0`$, as well as its first derivative, must be zero. As a function of $`p_t`$, the pions exhibit an almost linear dependence of $`v_2`$, whereas the protons + antiprotons exhibit a more quadratic behavior, clearly different from the pions. As expected in a hydrodynamic picture, the kaons lie between the pions and the protons + antiprotons. The 90% purity of the protons + antiprotons in the $`0.8p_t0.9`$ range leads to a maximum +5% systematic error on the $`v_2`$ value in this bin. The observed behavior may be the result of the interplay between the mean expansion velocity, the elliptic component of the expansion velocity, and the thermal velocity of the particles. A similar effect was discussed for the case of directed flow in voloshinv1 .

The differential elliptic flow, $`v_2`$($`p_t`$), is plotted for pions for three different centrality selections in the upper panel of Fig. 3, and for protons in the lower panel. The open triangles represent the most central 11% of the measured cross section. The open squares correspond to 11 – 45% of the measured cross section, and the open circles correspond to 45 – 85%. The uncertainties on the points are statistical only. The systematic uncertainty is smallest for the centrality region with the best reaction plane resolution and is estimated to be 20% for the most central bin, 8% for the mid-central bin, and 22% for the most peripheral bin. At a given $`p_t`$, the more peripheral collisions have the largest value of $`v_2`$($`p_t`$), and $`v_2`$($`p_t`$) decreases for more central collisions. For all three of the centrality ranges, in the measured $`p_t`$ range, the $`p_t`$ dependence of $`v_2`$ for pions is approximately linear.

We have fitted the data with a simple hydrodynamic-motivated model. This model is a generalization of the Blast Wave model from Siemens ; pasi2  assuming the flow field is perpendicular to the freeze-out hyper-surface.

$$v_2(p_t)=\frac{_0^{2\pi }𝑑\varphi _b\mathrm{cos}(2\varphi _b)I_2(\alpha _t)K_1(\beta _t)(1+2s_2\mathrm{cos}(2\varphi _b))}{_0^{2\pi }𝑑\varphi _bI_0(\alpha _t)K_1(\beta _t)(1+2s_2\mathrm{cos}(2\varphi _b))},$$
(1)
where $`I_0`$, $`I_2`$, and $`K_1`$ are modified Bessel functions, and where $`\alpha _t(\varphi _b)=(p_t/T_f)\mathrm{sinh}(\rho (\varphi _b))`$, and $`\beta _t(\varphi _b)=(m_t/T_f)\mathrm{cosh}(\rho (\varphi _b))`$. The basic assumptions of this model are boost-invariant longitudinal expansion and freeze-out at constant temperature $`T_f`$ on a thin shell, which expands with a transverse rapidity exhibiting a second harmonic azimuthal modulation given by $`\rho (\varphi _b)=\rho _0+\rho _a\mathrm{cos}(2\varphi _b)`$. In this equation, $`\varphi _b`$ is the azimuthal angle (measured with respect to the reaction plane) of the boost of the source element on the freeze-out hyper-surface pasi2 , and $`\rho _0`$ and $`\rho _a`$ are the mean transverse expansion rapidity ($`v_0=\mathrm{tanh}(\rho _0)`$) and the amplitude of its azimuthal variation, respectively. With $`s_2=0`$, our equation reduces to Eq.2 from pasi2 . In Fig. 2, the fit to the minimum-bias data with $`s_2=0`$ is shown as the dotted lines. The poor fit shows that the data can not be described under the assumption of a spatially isotropic freeze-out hyper-surface in the transverse plane. This led us to generalize Eq.2 from pasi2  to the case of a spatially anisotropic freeze-out hyper-surface, introducing one extra parameter, $`s_2`$, describing the variation in the azimuthal density of the source elements, $`2s_2\mathrm{cos}(2\varphi _b)`$. This additional parameter leads to a good description of the data, shown as the solid lines in Fig. 2. A positive value of the $`s_2`$ parameter would mean that there are more source elements moving in the direction of the reaction plane. The $`s_2`$ parameter does not distinguish between an in-plane or out-of-plane extended source, because the direction of the flow field is not an observable. However, azimuthally sensitive HBT measurements will be able to address this. In this model (Eq. 1), elliptic flow as a function of particle transverse momentum, $`v_2`$($`p_t`$), starts from zero and rises approximately quadratically with $`p_t`$ until the particle becomes relativistic, and then $`v_2`$($`p_t`$) continues to rise almost linearly. For the heavier particles, the linear rise is delayed. These mass-dependent effects are larger for lower temperatures ($`T_f`$) and larger transverse rapidities ($`\rho _0`$). The linear rise increases with both the amplitude of the azimuthal variation of the transverse expansion rapidity ($`\rho _a`$) and the elliptic deformation ($`s_2`$), but only $`\rho _a`$ produces a mass dependence. This arises from the dependence of the $`s_2`$ parameter effect on momentum not energy. This also results in a change in slope and a characteristic curvature in $`v_2`$($`p_t`$), since the flow associated with $`s_2`$ saturates quite early, when $`p_t/T1`$ pasi3 . The fact that the data can not be described with $`s_2=0`$ leads to the interpretation that the elliptic flow is not caused by an azimuthal velocity variation alone, but by the combination of a velocity difference and a spatially anisotropic freeze-out hyper-surface.

Table 1 lists the results of the two fits (without, and with $`s_2`$) for minimum-bias data. The parameters are correlated and therefore have large uncertainties. It is expected that better constraints on the parameters can be obtained in the future, when the fits to the single inclusive spectra for different mass particles are available. As the solid lines in Fig. 2 show, the data could be described by a reasonable set of parameters. For the different centralities, the $`T_f`$ and $`\rho _0`$ were kept fixed to the minimum-bias values. Both the obtained $`\rho _a`$ and $`s_2`$, not shown, decrease from peripheral to central collisions as expected from a simple geometrical picture.

A full hydrodynamic calculation, shown in Ref. pasi2 , describes the $`v_2`$($`p_t`$) for pions, kaons and protons + antiprotons for minimum-bias collisions equally well. The agreement is consistent with early local thermal equilibrium and the presence of an early stage pressure gradient. However, within the current statistical uncertainties in the data and different theoretical interpretations teaney , the measurement does not justify drawing inferences about the different EOS. Comparing the centrality dependence from our Fig. 3 with the same hydrodynamic calculation, shows a deviation which is most pronounced for the most peripheral data. This might indicate that the amount of particle rescattering in peripheral collisions is insufficient to justify the local thermal approximation implicit in hydrodynamical models.

We have made the first measurement of identified particle elliptic flow at RHIC. The measured elliptic flow as a function of $`p_t`$ and centrality differs significantly for particles of different masses. This mass dependence can be described with a simple hydrodynamic-motivated Blast Wave model. This model suggests that elliptic flow is generated by the combination of an azimuthal velocity variation and a spatially anisotropic freeze-out hyper-surface. The mass dependence of $`v_2`$($`p_t`$) is also in close agreement with full hydrodynamic model calculations, suggesting that in Au+Au collisions at $`\sqrt{s_{_{NN}}}=130`$ GeV, a system is created which, for central and mid-peripheral collisions, is consistent with early local thermal equilibrium followed by hydrodynamic expansion.

Acknowledgments: We wish to thank P. Huovinen for providing us with the hydrodynamical model calculation, the RHIC Operations Group at Brookhaven National Laboratory for their tremendous support and for providing collisions for the experiment. This work was supported by the Division of Nuclear Physics and the Division of High Energy Physics of the Office of Science of the U.S.Department of Energy, the United States National Science Foundation, the Bundesministerium für Bildung und Forschung of Germany, the Institut National de la Physique Nucleaire et de la Physique des Particules of France, the United Kingdom Engineering and Physical Sciences Research Council, Fundacao de Amparo a Pesquisa do Estado de Sao Paulo, Brazil, and the Russian Ministry of Science and Technology.


text:
# Membership of 23 stars towards the bulge globular clusters NGC 6528 and NGC 6553 Observations collected at the European Southern Observatory - ESO, Chile

## 1 Introduction

The study of the stellar populations in the Galactic bulge is very important to constrain possible models of galaxy formation. In particular, the determination of the metallicities and abundance ratios of bulge stars, either from the field or in clusters, provides key information to help decide among the possible scenarios for the history of chemical enrichment of the Galaxy.

There are few previous studies of radial velocities and metallicity estimations of bulge stars from low resolution spectra. In the analysis of 400 field bulge stars by Sadler et al. (1996), metallicities and \[Mg/Fe\] values were estimated. Minniti (1995a,b) studied the membership of stars towards 7 bulge globular clusters.

The best studied among the bulge globular clusters are NGC 6528 ($`\alpha _{2000}`$ = 18h04m49.6s, $`\delta _{2000}`$ = -30o03’20.8”) and NGC 6553 ($`\alpha _{2000}`$ = 18h09m15.7s, $`\delta _{2000}`$ = -25o54’27.9”). Ortolani et al. (1995) have shown that, besides being old, these clusters have luminosity functions which are very similar to that of Baade’s Window, which indicates that they belong to the same stellar population.

NGC 6528 is located in the Baade Window, at a distance d = 7.83 kpc from the Sun, and NGC 6553 is relatively close to the Sun, at a distance d = 5.1 kpc (Barbuy et al. 1998). As they are both located in crowded fields, the measurement of radial velocities of individual stars is of crucial importance for the determination of their membership in the clusters.

Both clusters are known to be metal-rich. However, there is no consensus in the literature regarding their detailed metal abundances. Recently, Barbuy et al. (1999) analysed high resolution spectra of two giant stars of NGC 6553. An Iron abundance of \[Fe/H\] = -0.55$`\pm `$0.2 and abundance ratios \[Na/Fe\] $``$ \[Al/Fe\] $``$ \[Ti/Fe\] $``$ +0.5, \[Mg/Fe\] $``$ \[Si/Fe\] $``$ \[Ca/Fe\] $``$ +0.3 were derived. These ratios imply an overall metallicity \[Z/Z\] $``$ -0.1. Cohen et al. (1999), analysing high resolution spectra of five red horizontal branch stars, obtained a mean metallicity \[Fe/H\] = -0.16 and an excess of the $`\alpha `$-element calcium to iron of about 0.3 dex, which imply an overall metallicity \[Z/Z\] $``$ +0.1. Metal abundances of these clusters are also discussed in Barbuy et al. (1999) and Barbuy (2000).

In view of the disagreement between previous determinations of \[Fe/H\], it is important that abundance estimations be extended to a larger number of stars of both clusters. In this paper, we determine radial velocities, effective temperatures, gravities and estimations of metallicities \[Fe/H\] based on low resolution spectra for 23 stars towards NGC 6553 and NGC 6528, and verify their membership in these clusters.

In Sect. 2 the observations are described. The radial velocities derived are presented in Sect. 3. In Sect. 4 the stellar parameters are derived, and synthetic spectra are compared to observations to estimate metallicities. Concluding remarks are given in Sect. 5.

## 2 Observations

Low resolution spectra of individual stars of NGC 6528 and NGC 6553 were obtained in 1992 August and 1994 June, at the 1.5m ESO telescope at ESO (La Silla). The Boller & Chivens spectrograph was employed. In 1992 August the Thompson CCD # 18 with 1024$`\times `$1024 pixels, with a pixel size of 19 $`\mu `$m was used. A resolution of $`\mathrm{\Delta }\lambda `$ $``$ 8 Å and a spectral coverage of $`\lambda \lambda `$ 4800-8800 Å were achieved. In 1994 June, the Ford Aerospace FA 2048 L, frontside illuminated, uncoated CCD detector (ESO # 24) with 2048$`\times `$2048 pixels and pixel size 15$`\times `$15 $`\mu `$m was used. The grating # 27 resulted in a spectral resolution $`\mathrm{\Delta }\lambda `$ $``$ 4 Å and a spectral coverage in the range $`\lambda \lambda `$ 4800-7550 Å.

The log of observations is provided in Table 1. The stars are identified according to the charts by Hartwick (1975) for NGC 6553 and van den Bergh & Younger (1979) for NGC 6528. Spectra of a given star were co-added by weighting their S/N ratios; the final S/N are indicated in Table 1.

In Figs. 1 and 2 are shown the $`V`$ vs. $`VI`$ Colour-Magnitude Diagrams of NGC 6528 and NGC 6553 using data obtained with the Hubble Space Telescope (Ortolani et al. 1995) where the sample stars are identified.

## 3 Radial velocities

The radial velocities were determined by means of three methods, as explained below and reported in Table 2. The observed radial velocities derived were transformed to heliocentric values using the observation dates given in Table 1.

(a) A Fourier cross-correlation was applied on the program spectra relative to selected template spectra. As templates, 12 G, K and M stars were selected from the Jacoby et al.’s (1984) library, which have approximately the same spectral resolution (4.5 Å) of the sample spectra. The templates adopted were the ones closest in spectral type to each of the program stars. The spectra of both sample and template stars were normalized and different regions in the spectra were defined in order to give highest peaks of cross-correlation for each considered template. The results obtained with this method are given in column 2 of Table 2. The r.m.s. of the values derived with each template spectrum is of the order of 15 km s-1. A systematic effect was identified for the coolest stars, since all of them appeared to show lower velocities when compared to the hotter stars of the same cluster and these values were not considered.

(b) Mean shifts between the observed wavelengths of identified absorption lines and laboratory wavelengths were measured (column 3 of Table 2). The r.m.s. of the values derived is of the order of 15 km s-1.

(c) the code HALO (Cayrel et al. 1991) derives radial velocities by comparing the observed spectrum to a grid of synthetic spectra, using a cross-correlation technique. The grid of synthetic spectra available (Barbuy et al. 2001) does not contain stars cooler than Teff $`<`$ 4000 K, and for this reason the errors should be higher for velocities of stars cooler than Teff $`<`$ 3700 K in which TiO bands are pronounced.

Histograms of radial velocities of individual stars (coolest stars excluded) corresponding to each method were built. Gaussian curves were fitted to each histogram, from which the radial velocity corresponding to each method was derived, as reported in Table 3 together with values from the literature. An example of this procedure is presented in Fig. 3 for the cross-correlation technique using IRAF. The final radial velocities adopted for the clusters correspond to the mean of the values derived from the three methods.

## 4 Stellar parameters

### 4.1 Temperatures

The effective temperatures were estimated from $`BV`$, $`VI`$, $`VK`$ and $`JK`$ colours, based on the colour vs. Teff calibrations by Bessell et al. (1998), which in turn are based on NMARCS models by Plez et al. (1992) and their grid extensions. These effective temperatures are listed in Table 5. V and I colours were obtained with the Hubble Space Telescope (Ortolani et al. 1995) and J and K colours were obtained with the detector IRAC2 at the 2.2m telescope of ESO (Guarnieri et al. 1998).

For NGC 6528 the colour excesses adopted were E(V$``$I) = 0.68 and E(B$``$V) = 0.52 (Barbuy et al. 1998). For NGC 6553 E(V$``$I) = 0.95 and E(B$``$V) = 0.7 were adopted (Guarnieri et al. 1998). The $`VK`$ and $`JK`$ colours were dereddened assuming E(V-K)/E(B-V) = 2.744 and E(J$``$K)/E(B$``$V)=0.527 (Rieke & Lebofsky 1985).

An independent method for the derivation of temperatures was based on calibrations of equivalent widths of TiO bands. The indices as defined in Table 4 were measured on the grid of synthetic spectra by Schiavon & Barbuy (1999) in the range of parameters 3000 $`T_{\mathrm{eff}}`$ 5000 K, $``$0.5 $``$ log g $``$ 2.5 and \[Fe/H\] = $``$0.3. These indices are shown in Fig. 4 for a resolution of $`\mathrm{\Delta }\lambda `$ = 8 Å. Polynomial curves of the form Teff = f(W(TiO)) were derived and applied to the indices measured in the sample stars. The TiO indices are strongly sensitive to temperature for $`T_{\mathrm{eff}}`$ 3800 K as illustrated in Fig. 5. For $`T_{\mathrm{eff}}`$ 4000 K a degeneracy appears due to the fact that TiO bands are not present at these higher temperatures.

A more general polynomial of the form log W(TiO) = (a + b log Teff \+ c log $`g`$ \+ d \[Fe/H\] \+ e (log Teff)2 \+ f \[Fe/H\]2 \+ g \[Fe/H\] log Teff), valid in the range 2500 $``$ Teff $``$ 5000 K, -0.5 $``$ log g $``$ 2.5 and -0.5 $``$ Fe/H $``$ 0 was derived. The coefficients of the formula above are shown in Table 6 for convolutions of $`\mathrm{\Delta }\lambda `$(FWHM) = 4 Å and 8 Å.

The temperatures obtained and final values adopted are reported in Table 5. Photometric temperatures were adopted for stars for which the TiO temperature TTiO $``$ 3800 K, whereas for stars with TTiO $`<`$ 3800 K the mean of TiO temperatures were adopted.

### 4.2 Gravities

Gravities were derived using the classical relation log g= 4.44 + 4log T/T \+ 0.4(Mbol-Mbol⊙) \+ log M/M, adopting T = 5770 K, M = 0.8 M and Mbol⊙ = 4.74 cf. Bessell et al. (1998). For deriving Mbol∗ we used the distance modulus adopting a total extinction AV = 2.43 for NGC 6553 and AV = 1.8 for NGC 6528 (Barbuy et al. 1998). The bolometric magnitude corrections were taken from Bessell et al. (1998).

The resulting Mbol∗ and gravities are given in Table 7. Taking into consideration the errors due to uncertainties in Teff and Mbol the final error in log g is estimated to be of $`\pm `$0.5 dex.

### 4.3 Metallicities

Spectrum synthesis calculations were used to fit the observed spectra. The calculations of synthetic spectra were carried out using the code described in Barbuy et al. (2000) where molecular lines of MgH A2$`\mathrm{\Pi }`$-X2$`\mathrm{\Sigma }`$, CH A2$`\mathrm{\Delta }`$-X2$`\mathrm{\Pi }`$, CN A2$`\mathrm{\Pi }`$-X2$`\mathrm{\Sigma }`$, C2 Swan A3$`\mathrm{\Pi }`$-X3$`\mathrm{\Pi }`$ and TiO $`\alpha `$ C3$`\mathrm{\Delta }`$-X3$`\mathrm{\Delta }`$, $`\gamma `$ A3$`\mathrm{\Phi }`$-X3$`\mathrm{\Delta }`$ and $`\gamma `$’ B3$`\mathrm{\Pi }`$-X3$`\mathrm{\Delta }`$ systems are taken into account.

For atomic lines the laboratory oscillator strengths by Fuhr et al. (1988), Martin et al. (1988), Wiese et al. (1969), and laboratory values compiled by McWilliam & Rich (1994) were adopted whenever available, otherwise they were taken from fits to the solar spectrum (see discussion in Barbuy et al. 1999).

ATLAS9 and NMARCS models were employed. A grid of models using the ATLAS9 code (Kurúcz 1993) was created adopting a mixing length parameter $`\alpha `$ = 0.5 (see Barbuy et al. 2001). NMARCS photospheric models for giants by Plez et al. (1992) and their unpublished extended grids were employed (see more details in Schiavon & Barbuy 1999).

The metallicities were obtained based on two methods, both using synthetic spectra:

(i) The observed spectra were compared to synthetic spectra in the range $`\lambda \lambda `$ 5000–7500 $`\mathrm{\AA }`$. The metallicities were estimated by interpolating between synthetic spectra of \[Fe/H\]= 0.0, –0.3, –0.5 and –0.6, in all cases assuming \[Mg/Fe\] = +0.3, and the temperatures and gravities determined in Sects. 4.1 and 4.2.

(ii) Comparisons with a grid of synthetic spectra in the wavelength region $`\lambda \lambda `$ 4600-5600 Å, using the differences method as described in Cayrel et al. (1991) and Barbuy et al. (2001), are carried out. In this method, the observed spectrum is divided by a reference synthetic spectrum. The resulting signal can be expressed as a linear combination of variations in temperature, gravity and metallicity. In conjunction with the grid of synthetic spectra, it is possible to establish the differences in Teff, log g and \[Fe/H\] between the program star and the reference synthetic spectrum through a perturbation method. The grid covers the range 4000 $``$ Teff $``$ 7000 K, 0.0 $``$ log g $``$ 5.0, -3.0 $``$ \[Fe/H\] $``$ +0.3, and \[Mg/Fe\] = 0.0 and +0.4. Fig. 6 shows the fit to NGC 6528 I-1.

In Table 7 are listed the temperatures, gravities, \[Fe/H\] and \[Mg/Fe\] obtained with methods (i) and (ii). Note that method (ii) tends to give lower metallicities relative to method (i). This may be due to limitations of the grid of synthetic spectra, which is being extended to cover wider ranges of parameters.

The stars NGC 6528 I-5, NGC 6553 III-2 and II-95 are probable non-members, given that their atmospheric parameters are incompatible with their location in the Colour-Magnitude Diagrams of the clusters. The star II-51 appears to be too cool (Table 7) with respect to its CMD locus (Fig. 2). However, considering that it could be a red variable, this star is tentatively classified as a possible member.

Mean metallicities of \[Fe/H\] = -0.5 for NGC 6528 and \[Fe/H\] = -0.7 for NGC 6553 were estimated from gaussian fits to the histograms of metallicities given in Table 7 (typical standard deviations are of 0.2 dex). These metallicities, together with the magnesium excess of \[Mg/Fe\] $``$ +0.3, and assuming that all other $`\alpha `$ elements show an excess of +0.3 dex relative to iron, result in overall metallicities of \[Z/Z\] = -0.25 and -0.45 for NGC 6528 and NGC 6553 respectively.

## 5 Conclusions

The study of individual stars in globular clusters along their evolutionary stages is of prime importance for an improved understanding of stellar evolution. Low resolution spectroscopy provides a means for the study of a large number of stars. In the present work we have measured radial velocities and estimated metallicities in 23 stars towards the globular clusters NGC 6528 and NGC 6553, which allows us to identify member stars. We also obtained their atmospheric properties to a first approximation. This is an important step before applying efforts to obtain high resolution spectroscopy with 8m class telescopes. The method presented here is also of interest for last generation multi-object instruments such as VLT-VIMOS.

The stars were analysed by comparisons between their observed spectra and a grid of synthetic spectra. TiO equivalent widths were used to estimate effective temperatures of stars cooler than $`T_{\mathrm{eff}}`$ $``$ 3800 K and a calibration of equivalent widths of TiO bands as a function of atmospheric parameters is presented.

Mean values of heliocentric radial velocities of vr = 234 km s-1 for NGC 6528 and vr = -1 km s-1 for NGC 6553 are derived.

Regarding membership, among the 23 stars observed we concluded that 4 of them are probable non-members. These are: NGC 6528 II-39, non-member due to a deviant radial velocity, and NGC 6528 I-5, NGC 6553 III-2 and II-95, non-members due to incompatibilities of atmospheric parameters vs. location in the Colour-Magnitude Diagrams.

NGC 6553 II-51 could be a non-member, or a red variable for which the spectrum was taken during a cool phase.

The basic stellar parameters derived show the interesting result that there is a trend for member giants of NGC 6528 to be more metal-poor than the two Horizontal Branch stars NGC 6528 I-23 and I-24, thus reproducing the discrepancy found between analysis of NGC 6553 giants by Barbuy et al. (1999) and Horizontal Branch stars by Cohen et al. (1999). Given the errors involved in the analysis of low resolution spectra, these results have to be checked with high resolution spectra, and further studies of this discrepancy will be possible only with a homogeneous analysis of stars ranging from the red giant branch to the HB, and also employing different sets of model atmospheres all along the evolutionary sequence.

In summary, we obtained for NGC 6528 and NGC 6553 metallicities of \[Fe/H\] = -0.5 $`\pm `$ 0.3 and \[Fe/H\] = -0.7 $`\pm `$ 0.3. Using \[Mg/Fe\] $``$ +0.3, and assuming that other $`\alpha `$ elements show the same excess of +0.3 dex relative to iron, the results are \[Z/Z\] = -0.25 and -0.45 for NGC 6528 and NGC 6553 respectively.

###### Acknowledgements.

We are grateful to A. Milone for having carried out part of the observations. We acknowledge partial financial support from CNPq and Fapesp. P. Coelho and T. Idiart acknowledge respectively the Fapesp Master fellowship no 98/07492-4, and Post-Doc fellowship no 97/13083-7. RPS acknowledges support provided by the National Science Foundation through grant GF-1002-99 and from the Association of Universities for Research in Astronomy, Inc., under NSF cooperative agreement AST 96-1361.


text:
# Discrimination of the Bell states of qudits by means of linear optics

## I Introduction

Quantum optics and particularly linear optical elements represent important tools for the experimental investigation of the basic features of quantum information transfer and processing and fundamentals of quantum theory. The special interest is devoted to the study of entangled states that are not only highly interesting by themselves but that find also the use in quantum teleportation, quantum dense coding, quantum cryptography, quantum computing, etc. Closely related to many of these issues is the so called Bell-state analysis, i.e., the discrimination of maximally entangled states completing the nonlocal orthonormal basis of two- or multi-partite quantum system.

The relative success concerning the discrimination of two of the four Bell states of two qubits with certainty by passive linear optical elements  challenges the question whether it is possible to do something similar also for the Bell states of two qudits if their dimension $`d>2`$. At the first glance the direct generalization of the approach given in Ref.  indicates that the limitations on the probability of success could be less restrictive for $`d>2`$ than for $`d=2`$. But when one starts to play with possible extensions of the original Innsbruck scheme  for qudits he quickly gets into troubles. The recent results of Calsamiglia  indicate that it is impossible without additional auxiliary photons. This was the motivation to subject the problem of “generalized” linear optics Bell-state measurement to more detailed analysis.

In the considered scheme, qudits are represented by $`d`$ modes of radiation (with the same frequencies and polarizations). The total number of photons in these $`d`$ modes is required to be one. The $`i`$-th “logical” (or computational) basis state corresponds to the situation when exactly one photon is present in the $`i`$-th mode. Such an implementation allows us to realize any unitary operation on a single qudit (up to a global phase) in a deterministic way by the means of passive linear optical elements .

In fact, we restrict our tools to beam splitters, phase shifters and delay lines. All of them may be electronically switched — the conditionally dynamics is allowed, i.e., some operations may depend on the result of the measurement on selected modes outgoing the “previous” operation. We consider ideal detectors that can distinguish the number of impinging photons (in a given mode). Even if such detectors are not available in practice yet our choice is justified as we seek for fundamental limitations of linear optical devices.

By the Bell states of two qudits we mean the maximally entangled states of the following form (see, e.g., Ref. )

$$|\psi _{mn}=\frac{1}{\sqrt{d}}\underset{j=0}{\overset{d1}{}}\mathrm{exp}\left[2\pi i\frac{jn}{d}\right]𝖺_j^{}𝖻_{jm}^{}|\text{vac},$$
(1)
where $`𝖺_k^{}`$ and $`𝖻_k^{}`$ are bosonic creation operators in corresponding modes of the first and the second qudit, respectively, $`|\text{vac}`$ denotes a vacuum state, and $`jm=(j+m)\text{ mod }d`$; $`m`$ and $`n`$ go from 0 to $`d1`$.

## II No-go theorem for the case without auxiliary photons

Calsamiglia has shown  that any linear optical device that does not use auxiliary photons in ancillary modes (however, that may include conditional dynamics) cannot unambiguously discriminate any state with Schmidt rank higher than two from any set of two-qudit states spanning the whole two-qudit Hilbert space (qudits are supposed to be represented by one-photon states over $`d`$ modes). This is true even if the probability of successful discrimination is allowed to be lower than one. If there is a nonzero probability of some detection event for some state from the set that has Schmidt rank higher than two then there is always at least one another state from the set that gives a nonzero probability for this detection event too.

It directly follows that there is no way how to distinguish any one of $`d^2`$ Bell states of two qudits (for $`d>2`$) without error by means of linear optics if no auxiliary photons are involved.

## III Impossibility to discriminate a Bell-state with certainty

In the paper of Carollo and Palma  it is shown under the same conditions as assumed here (i.e., linear elements, arbitrary number of auxiliary modes, conditional measurements, and photon number detectors are assumed) that any two $`L`$-photon states over $`M`$ modes, randomly chosen from a known set of $`K`$ states, are completely (i.e., with certainty) distinguishable in the presence of auxiliary photons only if they are completely distinguishable in the absence of auxiliary photons. It also means that if some two states are not completely distinguishable in the case when no auxiliary photons are involved, than they cannot be distinguished even if any finite number of auxiliary photons are employed.

As stated earlier it is not possible in any way to distinguish any one of $`d^2`$ Bell states from all of the remaining states by linear optical device with no additional photons if $`d>2`$. In other words, to each Bell state there is at least one another Bell state such that these two states cannot be distinguished from each other (providing $`d>2`$). According to the statement given in the previous paragraph this must stay valid even if auxiliary photons are allowed. So it is impossible to discriminate even one of $`d^2`$ Bell states with certainty (with 100 % probability of success) irrespective whether auxiliary photons are allowed or not.

## IV Efficient Bell-state analyzer with linear optics

Now we will show how to distinguish all the Bell states of two qudits with the probability of success arbitrarily close to one using only linear optical elements (and auxiliary photons).

Recently, a scheme for “non-deterministic” quantum computation with linear optics were proposed by Knill, Laflamme, and Milburn . It is based on a non-linear phase shift:

$$\alpha _0|0+\alpha _1|1+\alpha _2|2\alpha _0|0+\alpha _1|1\alpha _2|2,$$

where kets represent number states in a given mode. We will denote this operation “NS”. A simple device was designed that can perform this operation with probability 1/4. It consists of a couple of beam splitters and phase shifters. The effective non-linearity is provided by a measurement process. Two ancillary modes, with a single photon in one of them, and two photon-number detectors are necessary. The successful operations turns up when the first detector registers one and the second detector no photon.

Using two beam splitters and two operation NS one can built a conditional sign-flip gate (“C-SIGN”) as shown in Fig. 1. This gate inverts the sign of the state vector of two modes when there are exactly one photon in each of them. In the other cases (with at most one photon altogether) the operation does nothing:

$$\begin{array}{cccccccc}|1|1\hfill & & \hfill & |1|1,\hfill & |1|0\hfill & & & |1|0,\hfill \\ |0|1\hfill & & & |0|1,\hfill & |0|0\hfill & & & |0|0.\hfill \end{array}$$

The unitary matrix representing the transformation of creation operators on the beam splitter is supposed to be

$$\left(\begin{array}{cc}\mathrm{cos}\theta \hfill & \hfill \mathrm{sin}\theta \\ \mathrm{sin}\theta \hfill & \hfill \mathrm{cos}\theta \end{array}\right).$$

Particular angles $`\theta `$ are written inside the corresponding boxes in Fig. 1. Since this gate uses two NS operations its probability of success is 1/16.

However, in the mentioned papers the way is proposed how to increase the probability of the successful action of the C-SIGN arbitrarily close to one. This way is based on a “teleportation trick”. It requires preparation of relatively complicated ancillary entangled state of $`2n`$ photons in $`4n`$ modes:

$`|\varphi _n={\displaystyle \underset{j,k=0}{\overset{n}{}}}(1)^{(nj)(nk)}`$ $`|1^j|0^{nj}|0^j|1^{nj}`$ (2)
$`\times `$ $`|1^k|0^{nk}|0^k|1^{nk},`$ (3)
where $`|x^y=|x|x\mathrm{}|x`$, $`y`$-times. First $`n`$ modes are “put together” with one input mode and all these $`n+1`$ modes are subjected to $`n+1`$ point Fourier transform. Then the photon number measurement is performed on the transformed modes and according to the result one of the other $`n`$ modes is chosen as an output of the gate and its phase is modified in general. The same action is done with the next $`2n`$ modes and the other input mode; see Fig. 2. It can be shown that the total probability of the success of such a C-SIGN gate is

$$p=\left(\frac{n}{n+1}\right)^2.$$
(4)
The Fourier transform, selection of modes, and phase shifts can be implemented by linear optical elements in a deterministic way. The state (3) can be prepared by means of NS operations, beam splitters and phase shifters. The probability of successful preparation can be rather low. But we should stress that this concerns just the preparation of an ancilla. In principle, it can be being prepared in advance and one can try many times. For more details, including the estimation of the success probability of the preparation procedure and the number of necessary elements, see Refs. .

Having a C-SIGN gate one can realize a gate “C-SWAP” that conditionally swaps the two modes of radiation (provided there is at most one photon altogether in them). Its scheme is in Fig. 3. If there is one photon in the mode 1 the modes 2 and 3 are swapped. If no photon is present in the mode 1 they are not changed. The probability of the success of the C-SWAP is the same as for the C-SIGN. The gate C-SWAP is the key ingredient to construct a logical operation acting on two qudits $`x,y`$ of arbitrary dimension $`d`$ that we will call “C-SHIFT”:

$$xx,y(yx)\text{ mod }d.$$

For any $`x`$ this operation represents a cyclic permutation or “rotation” of the values of the second qudit. In total, $`d1`$ rotations are necessary — each of them corresponds to one possible value of the first qudit except the value zero that leads to identity transformation. Any such rotation can be implemented by at most $`d1`$ “transpositions”, i.e., C-SWAPs. Thus if the probability of success of the C-SWAP is $`p`$ then the total probability of the success of the C-SHIFT is at least $`p^{(d1)^2}`$. In particular cases it can be even better, e.q., for $`d=3`$ the probability of the success of the C-SHIFT is $`p^3`$. Corresponding “network” is in Fig. 4. If one uses C-SIGN gates described above with the probability of success given by Eq. (4) then the probability of the successful action of the C-SHIFT gate is

$$P\left(\frac{n}{n1}\right)^{2(d1)^2}.$$
(5)
Now we can build a Bell-state analyzer. To do it we need our C-SHIFT gate,

$$\text{C-SHIFT: }|x_1|y_2|x_1|(yx)\text{ mod }d_2,$$
(6)
followed by a generalized Hadamard transform acting on the first qudit,

$$\text{HAD: }|x_1\frac{1}{\sqrt{d}}\underset{k=0}{\overset{d1}{}}\mathrm{exp}\left[2\pi i\frac{kx}{d}\right]|k_1,$$
(7)
where $`|j_l`$ represents the $`j`$-th “logical” state of the $`l`$-th qudit. The complete setup is shown in Fig. 5. If the sequence of these two operations is applied on a Bell state

$$|\psi _{mn}=\frac{1}{\sqrt{d}}\underset{j=0}{\overset{d1}{}}\mathrm{exp}\left[2\pi i\frac{jn}{d}\right]|j_1|(j+m)\text{ mod }d_2$$
(8)
\[see also Eq. (1)\] then the output state of the two qudits is $`|n_1|m_2`$. Thus, if both the qudits are realized as described above the input Bell state can be determined by a simple photodetection.

The generalized Hadamard transform is the unitary transformation of one qudit. As mentioned earlier any such operation can be realized in a determistic way with passive linear optical elements (for our implementation of qudits) — simply by combination of beam splitters and phase shifters. The C-SHIFT gate described above is a non-deterministic gate. Its probability of successful action is given by Eq. (5). Clearly, that value determines also the probability of the successful discrimination of an unknown Bell state. Increasing $`n`$ this probability can be made arbitrarily close to one provided the dimension $`d`$ is fixed.

Let us note, that the described approach can be extended to the discrimination of “generalized Bell states” of $`N`$ qudits, i.e., the following maximally entangled states that complete an orthonormal basis in the Hilbert space of $`N`$ qudits:

$`|\psi _{k_1,k_2,\mathrm{},k_N}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{d}}}{\displaystyle \underset{j=0}{\overset{d1}{}}}\mathrm{exp}\left[2\pi i{\displaystyle \frac{jk_1}{d}}\right]`$ (9)
$`\times `$ $`|j_1{\displaystyle \underset{i=2}{\overset{N}{}}}|(j+k_i)\text{ mod }d_i,`$ (10)
where $`k_i=0,\mathrm{},d1`$. In this case the C-SHIFT operation must be applied $`N1`$ times: between the first and the $`N`$-th qudit, the first and the $`(N1)`$-th qudit, etc. Finally between the first and the second qudit. Then the generalized Hadamard transform must be performed on the first qudit. If there was a generalized Bell state (10) in the input then the output state reads

$$|k_1_1|k_2_2\mathrm{}|k_N_N.$$

## V Conclusions

We have interpreted some recent results concerning unambiguous state discrimination with linear optical elements from the point of view of Bell-state measurement in case of two qudits with dimension higher than two. This analysis leads to the conclusions that with no auxiliary photons it is impossible to discriminate such Bell states without errors and that it is impossible to discriminate such Bell states with certainty in any way by the means of linear optics.

On the other hand, we have shown by an explicit construction that it is possible, in principle, to build a linear optical Bell-state analyzer capable to discriminate all the states of the Bell basis of two (or even more) qudits with the probability of success arbitrarily close to one. This device is based on the generalization of the idea of linear optics quantum computation for qudits. The price for the high success probability is a complicated setup and a large number of required auxiliary photons in rather complex entangled states.

It seems that the key ingredient that is necessary for the increase of the probability of successful discrimination is the entanglement “added” through the ancilla. The methods of “linear optics quantum computation” enable us to prepare required entangled states from “separated” photons, in principle. But the probability of the success of such a preparation is very low. Besides this, the preparation of single photons represents itself a serious experimental problem as single-photon states are highly non-classical ones. By the way, recently it was shown for all pure input states and for the large class of mixed states that the beam splitter can serve for the preparation of an entangled state on its output only if the input state exhibits non-classical behavior .

## Acknowledgments

The author appreciates discussions with John Calsamiglia and Norbert Lütkenhaus. This research was supported under the European Union project EQUIP (contract IST-1999-11053) and the project LN00A015 of the Ministry of Education of the Czech Republic.