text:
# References

Influence of spin structures and nesting on Fermi surface and a pseudogap anisotropy in t-t’-U Hubbard model.

A.A.Ovchinnikov and M.Ya.Ovchinnikova.

Institute of Chemical Physics, RAS, Moscow, 119977, Kosygin str., 4.

## Abstract

Influence of two types of spin structures on the form of the Fermi surface (FS) and a photoemission intensity map is studied for $`tt^{}U`$ Hubbard model. Mean field calculations are done for the stripe phase and for the spiral spin structure. It is shown, that unlike a case of electron doping, the hole-doped models are unstable with respect to formation of such structures. The pseudogap anisotropies are different for $`h`$ and $`e`$ doping. In accordance with ARPES data for $`La_{2x}Sr_xCuO_4`$ the stripe phases are characterized by quasi-one-dimensional segments of FS in the vicinity of points $`M(\pm \pi ,0)`$ and by suppression of spectral weight in diagonal direction $`k_x=k_y`$. It is shown that spiral structures display the polarisation anisotropy: different segments of FS correspond to electrons with different spin polarisation.

PACS: 71.10.Fd, 74.20.Rp, 74.20.-z

The angular resolved photoemission spectroscopy (ARPES) is effective in the study of electronic structure of cuprates . It gives the image of the Fermi surface (FS) projection onto a $`CuO_2`$ plane. Conclusions of the early works (see , and references therein) were consistent with FS of the hole type centered at point $`Y(\pi ,\pi )`$ of two-dimensional Brillouin zone (BZ). Later other versions of the FS topology have been discussed. In particular, a presence of an electron type of the FS have been proposed with the center at $`\mathrm{\Gamma }(0,0)`$ in $`Bi_2Sr_2CaCu_2O_{8+\delta }`$ (BSCCO) . Revision of the problem  partly confirms the original version. At the same time for $`La_{2x}Sr_xCuO_2`$ (LSCO) it has been proved a transition from h-type of the FS to e-type during transition from underdoped (UD) to overdoped (OD) regions of the phase diagram . The ARPES data gave an evidence of the d-wave superconducting gap and of the opening of the pseudogap (PG) in UD BSCCO. Recently a bilayer splitting of bands  in BSCCO and the time reversal symmetry breaking in UD BSCCO  have been discovered.

A wide use of the photoemission intensity maps in the space $`k_x,k_y,\omega `$ poses a problem of extracting the FS from the ARPES data. One aspect of the problem, calculation of the matrix elements is discussed in . However the FS topology and intensities of various (main or shadow) segments of the FS depend also on the spin and charge structures of a ground state of a system.

The goal of the present work is a model study of the influence of periodic structures on the FS and the photoemission intensity based on the t-t’-U Hubbard model. Unlike the static structures in manganites, in cuprates one deals with rather dynamical structures with a time scale larger than time $`t>10^6÷10^9`$ sec in the $`\mu `$-SR experiments. So at small times or for the processes with energy resolution higher than $`\mathrm{}/t`$ the local spin structures may be considered in quasi-static limit. Then the question arises about a correspondence of the structures to the ARPES data. The experimental indications of the SDW and CDW structures in cuprates are the incommensurate peaks in the spin susceptibility at $`q=(\pi \pm \delta q,\pi ),(\pi ,\pi \pm \delta q)`$ in LSCO , linear structures along the $`CuO_2`$ bonds with period 4$`a`$ ($`a`$ is a lattice constant) observed in the tunnel spectra , periodic chess structure $`4a\times 4a`$ around the vortex in mixed state of BSCCO  etc.

Here we study the stripe and spiral structures in the simplest mean field (MF) approach, interpret the PG, discuss examples of structures with the different types of the PG anisotropy and to connect some properties of MF solutions with features observed in some cuprates. Earlier  the DDW (d- density wave) structure has been proposed as a possible hidden order parameter (OP) in cuprates. Here a search of possible OP is extended to stripe and spiral structures. Stabilization of such structures originates from the removing of a degeneracy of states in ”hot points” \- the van Hove singularities (VHS) or in parallel segments of the FS at the so-called nesting.

Consider a t-t’-U Hubbard model with the zero band energies

$`ϵ_k=2t(\mathrm{cos}k_x+\mathrm{cos}k_y)+4t^{}\mathrm{cos}k_x\mathrm{cos}k_y`$ (1)
Further we put $`t=1`$, so all energies are measured in units of $`t`$. In MF approximation this Hamiltonian is insufficient to describe the superconducting (SC) pairing. It may be supplemented by an empirical pairing interaction of electrons at neighboring centers deduced in more accurate approaches (for example, ). But here we study a normal state and therefore retain the original $`tt^{}U`$ model.

A periodic 2D structure is determined by two translation vectors

$`E_1=(E_{1x},E_{1y}),E_2=(E_{2x},E_{2y})`$ (2)
and the corresponding vectors of the inverse lattice $`B_1,B_2`$ which satisfy the relations $`E_iB_j=2\pi \delta _{ij}`$. (The components of $`E_i`$ $`B_i`$ are in units of the corresponding constants of the original and inverse lattices). The unit cell of the structure contains $`n_c`$ sites with coordinates $`j=(j_x,j_y)`$, so that each site $`n=n(L,j)=(n_1,n_2)=E_1L_1+E_2L_2+(j_x,j_y)`$ is determined by the integers $`L_1,L_2`$ determinig the position of the unit cell and by integers $`j=(j_x,j_y)`$, fixing the site inside the unit cell.

Let $`\stackrel{~}{k}`$ is the quasi-momentum inside the main BZ $`\stackrel{~}{G}`$ of a given periodic structure unlike the momentum $`k`$ from the BZ $`G`$ of the original lattice. The areas of $`\stackrel{~}{G}`$ and $`G`$ are restricted by conditions $`|\stackrel{~}{k}B_i|\pi `$ or $`|k_{x(y)}|\pi `$. The order parameters of the periodic structures are the mean charges and spin densities on sites of the unit cell

$`r_j={\displaystyle \frac{1}{N_L}}{\displaystyle \underset{L}{}}<r_{n(L,j)}>;S_{\mu j}={\displaystyle \frac{1}{N_L}}{\displaystyle \underset{L}{}}<S_{\mu ,n(L,j)}>`$ (3)
Here $`\mu `$ numerates the spin components, $`N_L=N/n_c`$ is a number of the unit cells, $`n_c`$ is a number of sites in it. In the MF approximation the mean energy (1) is $`\overline{H}=<T>+N_LU_j(r_j^2_\mu S_{j\mu }^2)`$ and the wave function is determined by the occupation of the one-electron eigenstates $`\chi _{k\lambda }^{}`$ of the linearized Hamiltonian

$$H_{lin}=T+N_L\underset{j}{}\{2Ur_j\widehat{r}_j2US_{\mu j}\widehat{S}_{\mu j}\}=\underset{\stackrel{~}{k}\stackrel{~}{G}}{}\widehat{h}_{\stackrel{~}{k}}.$$

The latter is divided into independent contributions from each reduced quasi-momentum $`\stackrel{~}{k}`$. Here $`\widehat{r}_j,\widehat{S}_{\mu j}`$ are the operators corresponding to averages (3). In the momentum representation the eigenstates of $`\widehat{h}_{\stackrel{~}{k}}`$ are expanded over a basis set of $`2n_c`$ Fermi operators

$`\chi _{\stackrel{~}{k}\lambda }^{}={\displaystyle \underset{m,\sigma }{}}c_{\stackrel{~}{k}+Bm,\sigma }^{}W_{m\sigma ,\lambda }(\stackrel{~}{k}),`$ (4)
where $`\lambda =1,\mathrm{},2n_c`$, $`Bm=B_1m_1+B_2m_2`$ and $`m=(m_1,m_2)`$ is such set of integers, that allows the vectors $`\stackrel{~}{k}+Bm`$ to share all phase space $`G`$.

The matrix of eigenvectors $`W_{m\sigma ,\lambda }`$ and the eigenvalues $`E_{\stackrel{~}{k},\lambda }`$ are determined by diagonalization of $`h_{\stackrel{~}{k}}`$ in the basis set $`\{c_{\stackrel{~}{k}+Bm,\sigma }^{}\}`$:

$`(h_{\stackrel{~}{k}})_{m\sigma ,m^{},\sigma ^{}}W_{m^{}\sigma ^{},\lambda }=W_{m\sigma ,\lambda }E_{\stackrel{~}{k},\lambda }`$ (5)
Here

$`(h_{\stackrel{~}{k}})_{m\sigma ,m^{},\sigma ^{}}=\delta _{mm^{}}\delta _{\sigma \sigma ^{}}ϵ_{\stackrel{~}{k}+Bm}+U{\displaystyle \underset{j}{}}\phi (j,m^{}m)[r_j\delta _{\sigma \sigma ^{}}S_{\mu j}(\sigma _\mu )_{\sigma \sigma ^{}}]`$ (6)
with $`\phi (j,m)=\mathrm{exp}[i(jj_0)Bm];Bm=B_1m_1+B_2m_2.`$ In turn, the order parameters (3) are expressed as

$`\{r_j,S_{\mu j}\}={\displaystyle \frac{1}{2N}}{\displaystyle \underset{\stackrel{~}{k}\stackrel{~}{G}}{}}{\displaystyle \underset{ms,m^{}s^{}}{}}\{\sigma _0,\sigma _\mu \}_{ss^{}}\phi (j,m^{}m)W_{ms,\lambda }^{}(\stackrel{~}{k})W_{m^{}s^{},\lambda }(\stackrel{~}{k})f(E_{\stackrel{~}{k}\lambda }\mu )`$ (7)
Here the Puali matrices $`\sigma _\mu ,\sigma _0`$ in (9) correspond to the components $`S_{\mu j},r_j`$ and $`f`$ is the Fermi function. Equations (9) determine the self-consistent MF solution.

The intensity of the photoemission of electron with a projection $`k`$ of the momentum on $`ab`$ plane and the energy $`E=h\nu \omega `$ is

$`I(k,\omega )=|M(k)|^2A(k\omega )f(\omega )R_{\omega k}.`$ (8)
It is determined by the matrix element $`M(k)`$, a spectral density $`A(k\omega )`$ and the Fermi function $`f`$. Usual convolution is done with the Gaussian function $`R_{\omega k}`$  with parameters, which characterize a finite resolution over $`k`$ and energy. The dependence of the matrix element $`M`$ on $`k`$ was studied in . Here for simplicity we take a constant value for it, since we study an effect of the structure on the spectral density $`A`$. In the one-electron approximation

$`A(k,\omega )={\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{~}{k}\stackrel{~}{G}}{}}{\displaystyle \underset{m,\sigma ,\lambda }{}}|W_{m\sigma ,\lambda }(\stackrel{~}{k})|^2\overline{\delta }(E_{\stackrel{~}{k}\lambda }\mu \omega )\delta _{k,\stackrel{~}{k}+Bm}`$ (9)
Here $`\lambda =1,\mathrm{},2n_c`$ and index $`m=(m_1,m_2)`$ numerates all independent Umklapp vectors $`Bm=B_1m_2+B_2m_2`$. A standard replacement of the $`\delta `$\- function in (9) by a function with finite width $`\mathrm{\Omega }`$ is implied. A map of $`I(k_x,k_y,\omega =0)`$ allows to visualize both the main and shadow segments of the FS. Though the band energies are periodic functions in $`k`$ space ( $`E_{\stackrel{~}{k}+Bm,\lambda }=E_{\stackrel{~}{k},\lambda }`$ for any $`m=(m_1,m_2)`$ ), the intensity (8) does not posses such periodicity. Therefore even for the constant matrix element $`M`$ in (8) various segments of the FS display themselves with a different intensity due to the Umklapp processes in the case of SDW or CDW structure.

The intensity map allows easily to repeat the known results  about changing of the FS topology with increase of the doping and opening the pseudogap in UD region for the homogeneous AF solutions. The energies of the lower and upper Hubbard bands are

$`ϵ_{\stackrel{~}{k}\lambda }=4t^{}\mathrm{cos}k_x\mathrm{cos}k_y\pm \sqrt{U^2d_0^2+4t^2(\mathrm{cos}k_x+\mathrm{cos}k_y)^2}+const`$ (10)
The VHS in the DOS of the lower band correspond to the ”hot points” $`M=(\pm \pi ,0),(0,\pm \pi )`$. A form of the FS critically depends on a sign of $`t^{}`$. At $`t^{}=0`$ energy $`E_\lambda (k)`$ is constant for all $`k`$ along the magnetic Brillouin zone (MBZ) boundary. At $`t^{}>0`$ the energies at the points M are lower than at the diagonal points. This leads to the formation of the known hole pockets around points $`(\pi /2,\pi /2)`$ or the electron pockets around $`(\pi ,0)`$, $`(0,\pi )`$ in UD case of the h- or e- doped models . In the same UD region the PG is open at $`x,y`$ directions or at diagonal directions for h- or e- doped models. Such a topology is determined by a profile of the energies $`E_\lambda (k)`$ of lower or upper Hubbard bands as a function of $`k`$ on the MBZ boundary. But this profiles are greatly influenced also by a formation of the spin or charge structures in system.

The MF calculations allow to test a stability or an instability of the homogeneous AF states with respect to the formation of the periodic spin or charge structures. The Fig.1 presents the mean energy $`\overline{H}=<H>`$ as a function of doping $`\delta =|1n|`$ for the normal state of the hole- and electron- doped models for a set of structures. In addition to homogeneous paramagnet (PM) and AF solutions the next structures are considered:

1) The stripe structure which consists of the antiphase AF stripes parallel to the y-axis. A fixed width 4a of the stripe (8 sites in a unit cell of the structure) is taken and the domain walls are chosen as centered on bonds or on the sites. Both structures have a similar form of the FS, but that with the bond-centered domain walls occurs to have some lower energy and only it is discussed below.

2) The spiral spin structures

$`<S_n>=d_0[𝐞_x\mathrm{cos}Qn+𝐞_y\mathrm{sin}Qn]`$ (11)
with spirality vector $`Q=Q_x=\pi (\eta ,1)`$ or $`Q=Q_{xy}=\pi (\eta ,\eta )`$, with $`\eta 0.750.8`$.

3) The checkered structure with antiphase AF square domain $`4a\times 4a`$. It have a higher energy than the above two structures, and unrealistic shape of the FS, so the results are not presented here.

4) The proposed in  states of the orbital antiferromagnet with staggered currents on a plackets. In MF approximation they do not appear for $`tt^{}U`$ Hubbard model without additional interactions. At hole doping and $`U/t5`$ there exist the MF solutions with staggered spin currents and rotations on $`\pi /2`$ of local spins of the neighbor sites. But their energies are higher than those of stripe and spiral structures and corresponding results are not presented here.

Most of calculations were carried out for structures with fixed commensurate period $`8a`$, though the optimal size of the AF domains or spiral vector depend on doping. For spiral states the dependence of $`Q(\delta )`$ has been calculated many times. Here we are interested only in the specific features of the FS and the pseudogap anisotropy.

Fig. 2 shows that at the h-doping the homogeneous AF solutions are unstable with respect to formation of the stripe and spiral structures. Formation of these structures extend the doping range where the local spin retains a nonzero value (a disappearing of the $`<S_{\mu j}>`$ corresponds to matching of the energy of some structure and that of the paramagnet state). At the same time at e-doping the MF energies of the above structures appear to be higher than the energy of the homogeneous AF solution. Large stability of the AF state in the e-doped models corresponds to the more wide region of the AF order in cuprates of the n-type than in ones of the p-type. The result about the absence of the stripe and spiral structures in the n-type models are consistent also with a commensurate peak at $`Q(\pi ,\pi )`$ in the spin susceptibility observed in NCCO, PCCO (unlike the incommensurate peaks in the hole doped cuprates).

Fig.2a presents a map of the intensity $`I(k,\omega =0)`$, Eqs.(8,9), for the stripe structure in the model with $`U/t=6`$, $`t^{}/t=0.1`$ at $`n=0.8`$. In case of AF MF solution the chosen doping is close to the ”optimal” doping at which the Fermi level passes through the VHS in the DOS of the lower Hubbard band. The stripes (which are parallel to an y-axis) split the VHS. The ”hot points” $`M_x=(\pi ,0)`$ and $`M_y=(0,\pi )`$ become non-equivalent points and the quasi-1D asymmetric horizontal segments of the FS are formed. Fig.2b presents the intensity map symmetrized over the structures with the x- or y- stripe orientation. The FS form differs strongly from that for the homogeneous AF state of the model with $`t^{}>0`$. The main difference is the absence of the FS in the diagonal direction and near the $`M_y`$. This gives an evidence that the stripes open the pseudogap in these regions of $`k`$. The corresponding band energies $`E_\lambda (k)`$ descend below the Fermi level and the work function $`\mathrm{\Delta }_{PG}=\mu E_\lambda >0`$ for the electron to be removed from such $`k`$ is observed as the PG in ARPES data. The revealed anisotropy of the FS and pseudogap may change our conclusions about possible symmetry of the SC order in the striped state. The retaining of the quasi-1D segments of the FS only suggests the 1D character of conductivity along the stripe direction (here y-axis). With the increase of $`t^{}`$ up to the value 0.3 there appear the small additional hole pockets at around $`k_S=(\pi /2,\pi /2)`$ in the intensity map.

The anisotropy of the FS and the PG suggests the revision of the possible symmetries of SC order coexisting with the stripes. If the structure is symmetric with respect to replacement $`xy`$, then the d-wave SC order $`<c_k^{}c_k^{}>\mathrm{cos}k_x\mathrm{cos}k_y`$ is expected. It provides the orthogonality of the pair function to the on-site pair function $`c_n^{}c_n^{}`$ which is suppressed by the on-site repulsion U. Taking into account a new FS view (see Fib.2a) one may suppose the extended s-wave SC pair function $`<c_k^{}c_{k,}^{}>\mathrm{cos}k_x+\mathrm{cos}k_y`$. The latter might be orthogonal to the on-site pair function due to the node lines $`k_x\pm k_y=\pm \pi `$. Unlike the gain of the node lines in the angular dependence of the d-SC gap, one can suppose a gain of the node in an analogue of a ”radial” part of the pair function. Verification of such hypothesis requires further calculations.

Fig.3a presents the one-electron energies $`E_\lambda (k)`$ (eigenvalues of MF problem) as functions of quasi-momentum varying along the contour

$$Y(\pi ,\pi )M_y(0,\pi )Y(\pi ,\pi )M_x(\pi ,0)Y(\pi ,\pi ).$$

As expected, band energies $`E_\lambda (k)`$ are periodic functions with the period $`\pi /4`$ on the first horizontal segment of contour $`YM_yY`$. However, the intensity map selects only the energy levels which are non-shadow at given $`k`$. Fig.3b presents such map of $`I(k,\omega )`$ on the plane $`k,\omega `$ for $`k`$ varying along the same contour. Intersections of band energies with the Fermi level in the vicinity of $`M_x`$ correspond to quasi-1D FS of Fig.2a.

The above pictures may be explained by an action of the spin-dependent average field. The main harmonic of this field is $`F(n)\mathrm{cos}\{Q_\eta n+\phi _0\}`$ with $`Q_\eta =(\eta \pi ,\pi )`$ (here $`\eta =0.75`$). This field shifts upstairs the zero level $`ϵ(0,\pi )`$ at $`M_y`$ pushing it from the lower levels $`ϵ_k`$ at $`k=(\pm \eta \pi ,0)`$ near $`M_x`$. The same field descends downstairs (lower the Fermi level) the zero level $`ϵ(\pi ,0)`$ at $`M_x`$ pushing it from the higher zero levels $`ϵ(\pm 0.25\pi ,\pi )`$ near $`M_y`$.

The ARPES data for UD LSCO  are in qualitative agreement with the properties of the model calculations. It was shown  a presence of two different segments of the FS - near the points $`M`$ and at diagonal directions, systematic suppression of spectral weight near $`(\pi /2,\pi /2)`$ in comparison with BSCCO or OD LSCO, the straight segments of the FS near $`(\pi ,0)`$ with the width $`\pi /2`$. All these observed features have been interpreted as an evidence of nonhomogeneous stripe structures in UD LSCO  or the combine action of the stripe structure and of the order-disorder . Our calculations confirm this interpretation. Another evidence of the stripe structure is given by the observation of the ”incommensurate” peaks in neutron scattering at $`k=(\pi \pm \delta ,\pi )`$, $`(\pi ,\pi \pm \delta )`$ in LSCO .

Another type of the spin and angular anisotropy of the FS manifests itself in spiral spin structures (11) which possess the polarisational anisotropy of the FS. Different segments of the FS correspond to the electrons with different preferential spin polarisation. Really, the average field from spiral spin structure mixes one-electron states $`\{c_{k,}^{},c_{k+Q,}^{}\}`$ provides a splitting of the VHS’s. But unlike the stripe structures, the occupations $`n_{k,\sigma }=<c_{k\sigma }^{}c_{k\sigma }>`$ and intensities $`I_\sigma (k\omega )`$ of the photoemission of the electrons with fixed polarisation $`\sigma `$ depends on this polarisation $`\sigma `$. Values $`I_\sigma (k\omega )`$ are determined by Eqs.(8,9), but without summing over $`\sigma `$ in right part of (9).

Fig.4a presents a map of the intensity $`I_{\sigma =}(k,\omega =0)`$ for the spiral state with $`Q=(0.8\pi ,\pi )`$ and for the spin polarisation $`\sigma =`$ on the axis $`z`$ normal to the plane of the spin rotation of the spiral structure (11). Thus a photoemission of electrons with the up spin polarization ($`\sigma =`$) corresponds to the strongly anisotropic FS. For the opposite spin polarisation the FS coincides with the FS of Fig.3a, reflected in the plane $`x=0`$. Fig.3b presents the FS, symmetrized over the spins and over two types of structures with $`Q=(0.8\pi ,\pi )`$ $`Q=(\pi ,0.8\pi )`$. In the vicinity of the point $`M`$ the FS have both the intersection with the line $`MY`$ typical for the FS of the h-type and the same with the line $`M\mathrm{\Gamma }`$ typical for the e-type of the FS. Such double intersections have been probably observed in ARPES data of BSCCO . Direct comparison with experiment would require to take into account the bilayer structure of BSCCO and the corresponding band splitting.

Polarisational anisotropy of the FS directly reflects the existence of the spin currents $`J_{}=J_{}`$ in the spiral state. According to  it might be a reason of the time reversal symmetry breaking observed in the UD BSCCO in dichroism of the ARPES signal . In the alternative hypothesis  this effect is explained by the specific alignment of the circular micro currents.

A direct observation of the polarisational anisotropy of the FS needs the selective over the spin polarisation measurements of the photoemission intensity. Such selective measurements of the total photoemission have been achieved in the so called ”spin-orbit” photoemission . In principle, similar selectivity is possible in ARPES also. It is important to continue the study of the time reversal symmetry breaking, in particular, to clarify, whether this effect and accompanying currents have the surface or bulk character. Yet we cannot answer the question: if the ground state of the BSCCO might be presented as a set of quasi-static domains with the spiral structure and corresponding systems of spin currents.

In conclusion, the MF treatment of the normal state of the $`tt^{}U`$ Hubbard model shows that the FS topology and the PG anisotropy depends on sign of $`t^{}`$, on the type (e- or h- ) of the doping and of the spin structure. For the hole-doped models the homogeneous AF MF state occurs to be unstable with respect to formation of the stripe phase and the spiral spin structure. However, in electron-doped models the lowest energy refers to the homogeneous AF solution. This corresponds to a wide doping range of an existence of the local AF order in the n-type cuprates NCCO, PCCO and to commensurate peak at $`Q=(\pi ,\pi )`$ observed in neutron scattering in these cuprates. In accordance with the ARPES data for the LSCO, the stripe phase is characterized by the quasi-1D behavior of the FS in the vicinity of the ”hot points” $`(\pm \pi ,0)`$, by the pseudogap opening and by suppression of the spectral weight in the diagonal direction and in the direction parallel to stripes. Such anisotropy of FS and PG does not conform with the d- symmetry of the SC order. For the spiral spin structure the polarisation anisotropy of FS is revealed which means that different segments of the FS correspond to different spin electron polarisations.

Work is supported by the Russian Foundation of Basic Research under Grant No.03-03-32141. Authors are greatfull to V.Ya.Krivnov for stimulating discussions.

㭪 .

Fig.1.

In the top fig. the doping dependencies of the mean energy (per one site) for the hole doped $`tt^{}U`$ model with parameters $`U=4.0`$, $`t^{}=0.3`$ for different spin structures are shown: for the paramagnet state (curve PM), for the homogeneous AF state (curve AF), for the stripe phase with a period of 8a along the x-axis (dashed line 1), for the spiral states with vectors $`Q=\pi (\eta ,1)`$ (curve 2) and $`Q=\pi (\eta ,\eta )`$ (curve 3) at $`\eta =0.8`$. The bottom graphics- the energies of the same structures (with the same notations) for the electron-doped model. For convenience the common function $`F(\delta )=U(n^21)`$ is subtructed from all energies.

Fig.2.

a. A map of the photoemission intensity $`I(k,\omega =0)`$ for the antiphase AF stripes parallel to the $`y`$ axis, having a width 4a and the bond-centered domain walls. b. The same, but averaged over two stripe phases of $`x`$\- and $`y`$\- orientation. Parameters of model are $`U=6`$, $`t^{}=0.1`$.

Fig.3.

a. The eigenenergies $`E_\lambda (k)`$ of the MF problem as functions of the quasi-momentum running along the contour $`Y(\pi ,\pi )`$-$`M_y(0,\pi )`$-$`Y(\pi ,\pi )`$-$`M_x(\pi ,0)`$-$`Y(\pi ,\pi )`$ for the stripe structure with period 8a along the x-axis. b. The map of the intensity $`I(k,\omega )`$ for $`k`$, running along the same contour. The map displays the same levels $`E_\lambda (k)`$, but with a spectral weight, determined by the structure of band states. Parameters of the model are the same, as in Fig.2.

Fig.4.

a. The map of the photoemission intensity $`I_{}(k,\omega =0)`$ for the electrons with a spin polarisatin $`\sigma =`$ for the spiral structure with $`Q=\pi (0.8,1)`$. b. The same but averaged over two spin polarisations $`\sigma =,`$ and over two structures with $`Q=(0.8\pi ,\pi )`$ and $`(\pi ,0.8\pi )`$. The model parameters are $`U=6`$, t’=0.1.


format: text

text:
# Description of the particle ratios and transverse-momentum spectra for various centralities at RHIC in a single-freeze-out model

## 1 Introduction

Ultra-relativistic heavy-ion collisions offer a unique possibility to produce hot and dense hadronic matter in laboratory conditions. With the new results presently coming from RHIC  we have entered an exciting era when various model predictions and theoretical ideas may be confronted directly with the large amount of diverse and accurate data. One of the key issues in the field is the question of the degree of thermalization of the hadronic matter produced in heavy-ion collisions. From the outset of the heavy-ion physics the thermal models have been useful and successful to a large degree in the description of data  (for a recent review see ).

In this paper we analyze the hadron production at mid-rapidity measured at RHIC at the full beam energy of $`\sqrt{s_{NN}}`$=200 GeV within the single-freeze-out model of Ref. . The new element of this work is the application of the approach to data at various values of the centrality parameter, $`c`$ (or the impact parameter, $`b`$), which allows us for a closer look at the geometry and flow of the hadronic system formed in non-central collisions. The obtained systematics for the geometric parameters as functions of $`c`$ is very interesting, as it conforms to the wounded-nucleon scaling . We find that the invariant time at freezeout, $`\tau `$, and the transverse size, $`\rho _{\mathrm{max}}`$, behave in a very systematic way for not-too-large values of $`c`$. In particular $`\rho _{\mathrm{max}}(1c)`$, and $`\tau /\rho _{\mathrm{max}}\mathrm{const}`$.

The single-freeze-out model was formulated in Ref.  where it was successfully used to describe the transverse-momentum spectra of hadrons measured at the then-available lower beam energy of $`\sqrt{s_{NN}}`$=130 GeV. The model combines a typical thermal approach, used to study the ratios of hadron abundances  with a hydrodynamic expansion including in a natural way the longitudinal and transverse flow. A characteristic and important feature of the model is the complete treatment of the resonances in both the calculation of the ratios of hadron yields and the analysis of the spectra. The model and the data at $`\sqrt{s_{NN}}`$=130 GeV, including strange-quark particles such as $`\varphi `$, $`\mathrm{\Lambda }`$, $`\mathrm{\Xi }`$, $`\mathrm{\Omega }`$, as well as $`K^{}`$, have been found to be in a surprising agreement . A comparison of the model predictions with the data for particle abundances and transverse-momentum spectra, collected at the maximum RHIC energy of 200 GeV and presented in detail in this paper, bring further evidence for rapid thermalization, or a statistical nature, of the hadronic matter produced at mid-rapidity in heavy-ion collisions. Moreover, with our simple approach we form a basis for the explanation of most of the soft features of the hadron production observed at RHIC.

The paper is organized as follows: In the next Section we outline the main assumptions of the model, bringing up for completeness some of the basic material of Refs. . In Sect. 3 we study the ratios of hadron yields at $`\sqrt{s_{NN}}=200`$ GeV and find the thermodynamic parameters characterizing the freeze-out. In Sect. 4 we calculate the transverse-momentum spectra. The calculations are done for different centrality windows, such as defined by the experiment. Naturally, since the geometry of the hot system formed in the collision depends on the centrality, we find a dependence of the geometric/expansion parameters on $`c`$. Interestingly, this dependence reflects quite accurately, for not too large $`c`$, the scaling of hadron production with the number of wounded nucleons.

## 2 The single-freeze-out model

In our previous work  we have shown that with the complete treatment of hadronic resonances, the distinction between the traditionally considered two freeze-outs, the chemical and the thermal (or kinetic) one (see the following discussion), is not necessary. At least for the RHIC energies one can achieve a very good explanation of the soft part of the hadronic data by assuming a single freeze-out which takes place at the universal temperature

$$T_{\mathrm{chem}}=T_{\mathrm{kin}}T.$$
(1)
Clearly, a theoretical description of the whole duration of the collision process of the two ultra-relativistic heavy ions is difficult, since different degrees of freedom are important at various stages of the evolution. Consequently, one is forced to describe each stage within a different theoretical framework. In our approach we use only hadronic degrees of freedom and thus concentrate on the latest stages of the evolution of the system, namely when the hadrons, at some instant, cease to interact and stream freely to the detectors. This moment is called the freeze-out. Admittedly, this is a far-reaching simplification. Generally speaking, the freeze-out in itself may be a complicated process, involving duration in time, and a hierarchy where different kinds of particles and different reactions switch off at different times. In general, by kinetic arguments, one expects that reactions with lower cross sections switch off at higher densities/temperatures, while those with larger cross sections last longer. Since, in most cases, the elastic cross sections are larger than the inelastic ones (a counter example, though, is the $`p\overline{p}`$ interaction), one expects that the inelastic (or chemical) freeze-out occurs earlier than the thermal freeze-out. On these grounds one may also argue that strange (or charmed ) particles decouple earlier than other hadrons. With such a picture in mind one may “separate” the freeze-out process into a series of more specific freeze-outs associated to particular reaction channels. In fact, the works of Heinz and collaborators  introduced the concept of the chemical and the thermal (or kinetic) freeze-outs, with $`T_{\mathrm{chem}}>T_{\mathrm{kin}}`$. At the chemical freeze-out the inelastic interactions between the hadrons cease, and the chemical composition of the system is fixed. Later, only the elastic interactions are effective, leading to further cooling of the system. At the thermal freeze-out even the elastic interactions become ineffective and the hadrons are completely decoupled from each other.

One should bare in mind, however, that the scales in the system responsible for the thermalization and expansion are of similar order. To be more precise, for expanding systems the criterion for the thermal freeze-out to occur is

$$\xi =\tau _{\mathrm{exp}}/\tau _{coll}1,$$
(2)
with

$$\tau _{\mathrm{exp}}=1/(_\mu u^\mu ).$$
(3)
With our parameterization for the expansion profile we obtain $`\tau _{\mathrm{exp}}10\mathrm{GeV}^1`$, while Ref.  gives $`\tau _{coll}10\mathrm{GeV}^1`$ at $`T165\mathrm{MeV}`$. Thus the discussion on the validity of our approach (early thermal freeze-out) as opposed to approaches having later thermal freeze-out becomes a discussion concerning the precise value of the factor in the criterion of Eq. (2). It is possible that already at the moment of hadronization the system is so dilute and possesses such a strong flow that not only inelastic, but even the elastic processes can be neglected. Thus our main assumption of Eq. (1) may be, in the least, a good starting point. At any case, assumptions of a phenomenological model like ours are verified a posteriori by the experiment. As we have shown in Refs.  and will also demonstrate in this work for the full RHIC energy of $`\sqrt{s_{NN}}=200`$ GeV, we are in a rather comfortable situation. The single-freeze-out hypothesis, when combined with expansion, offers a very economic description of the data, with very few parameters (two universal thermal parameters and two geometric parameters for each centrality bin).

An experimental argument in favor of a short time between the chemical and kinetic freeze-outs has been provided by the observation of resonances, such as the $`K^{}`$ or $`\rho `$, by the STAR Collaboration , in the invariant-mass spectra. In the thermal approach the abundance of these resonances is compatible with a high freeze-out temperature, around 160 MeV. On the other hand, if the elastic scattering processes ($`K\pi `$ for the $`K^{}`$ or $`\pi \pi `$ for the $`\rho `$) took place for a much longer time, such that the system would cool off significantly, the abundance of the resonances would be much lower. We should point out that the scenario of a single freeze-out, put forward in Ref. , is compatible with the sudden hadronization scenario of Ref. .

Thus, our starting point is the hypothesis (1). Similarly to other thermal models of hadron production, we fix the value of $`T`$, and also the value of the baryon chemical potential $`\mu _B`$, by fitting the ratios of hadron abundances. The results of such an analysis will be presented in the next Section.

The assumption of the single freezeout allows us to calculate uniformly, i.e. within the same model, both the ratios of hadron abundances and the transverse-momentum spectra. In order to do this, however, we need to make a choice of the freeze-out hypersurface (i.e. a three-dimensional volume in the four-dimensional space-time) and of the four-velocity field of expansion at the freeze-out. A priori, since we do not consider dynamics of the earlier stages here and only make an educated guess, many choices are possible and for each particular case one may check the validity of the proposed form of the hypersurface and flow by confronting the output of the model calculations to the data. Our choice has been made in the spirit of the hydrodynamic calculations of Ref. , and is defined by the condition

$$\tau =\sqrt{t^2r_x^2r_y^2r_z^2}=\mathrm{const}.$$
(4)
To make the transverse size of the fire-cylinder,

$$\rho =\sqrt{r_x^2+r_y^2},$$
(5)
finite, we impose the condition

$$\rho <\rho _{\mathrm{max}}.$$
(6)
In addition, we assume that the four-velocity of the hydrodynamic expansion at freeze-out is proportional to the coordinate (Hubble-like expansion),

$$u^\mu =\frac{x^\mu }{\tau }=\frac{t}{\tau }(1,\frac{r_x}{t},\frac{r_y}{t},\frac{r_z}{t}).$$
(7)
Such a form of the flow at freeze-out, as well as the fact that $`t`$ and $`r_z`$ coordinates are not limited and appear in the boost-invariant combination in Eq. (4), imply that our model is boost-invariant. In view of the recent data delivered by BRAHMS , this assumption turns out to be justified for the description of particle production in the rapidity range $`1<y<1`$, where the variation of the observed particle multiplicities is moderate.

We note that the freeze-out hypersurface and the flow are jointly controlled by the parameters $`\tau `$ and $`\rho _{\mathrm{max}}`$. Hence, our model has altogether only four parameters (dependent on the colliding energy and the centrality parameter), namely, the two thermodynamic parameters, $`T`$ and $`\mu _B`$, and the two expansion (geometric) parameters, $`\tau `$ and $`\rho _{\mathrm{max}}`$. The thermodynamic parameters turn out to depend on the centrality bin very weakly, reflecting the very weak observed dependence of the particle ratios. On the other hand, the geometric parameters do depend on $`c`$ and their dependence, obtained from independent fitting, assumes a very intuitive form (cf. Sect. 6).

Other choices for the freezeout hypersurface and expansion have been investigated in the literature , with a recent study made in Ref. , where a class of parameterizations has been tested with a limited number of resonances included.

All resonances from the Particle Data Tables  have been included in our study. The method and necessary formulas that describe the decays in cascades can be found in Ref. . We stress the important role of including the high-lying states, which although thermally suppressed, are increasingly numerous according to the Hagedorn hypothesis . Moreover, the inclusion of resonances is crucial for the description of the $`p_{}`$-spectra, as it increases the logarithmic slope as if lowering the temperature . This observation works in favor of the hypothesis (1).

## 3 Particle ratios

For boost-invariant systems the ratios of hadron abundances at mid-rapidity are equal to the ratios of hadron densities (for a more detailed discussion of this point, which includes the effects of the resonance decays, see Ref. ). Using this property, we can determine thermodynamic parameters of our model directly from the study of the ratios of hadron abundances. We stress that we have only two independent thermodynamic parameters in our model: the temperature, $`T`$, and the baryon chemical potential, $`\mu _B`$. The other two chemical potentials ($`\mu _S,\mu _I`$) are fixed by the conditions that the initial strangeness of the system is zero and that the ratio of the baryon number to the electric charge is the same as in the colliding nuclei. This assumption could be lifted at the expense of carrying an extra parameter. However, in practice this is irrelevant, since the value of $`\mu _I`$ turns out to be very small and the effects of the isospin violation can be neglected.

In Table 1 we show the results of our fit to the RHIC data collected at $`\sqrt{s_{NN}}`$= 200 GeV. We include here only the stable hadrons (with respect to strong interactions) and use the most central data. The values of $`T`$ and $`\mu _B`$ shown in Table 1 have been fitted with the $`\chi ^2`$ method to the experimental data listed in the third column. We stress a very good quality of the fit, reflected by the low value of $`\chi ^2`$ per degree of freedom. An interesting feature of our fit is the stable value of the temperature, when compared to other fits done at lower energies. The optimum value of $`T`$ is practically the same as that found in the analysis of the RHIC data at $`\sqrt{s_{NN}}`$= 130 GeV, and (within errors) it is the same as the temperature found at the top-energy collisions at CERN SPS .333The similarities of the SPS and RHIC spectra are discussed in Ref. . On the other hand we observe a clear drop of the baryon chemical potential, from 41 MeV found at $`\sqrt{s_{NN}}`$= 130 GeV down to 29 MeV found in the present calculation. A smaller value of the baryon chemical potential is a simple consequence of the observed increase of the $`\overline{p}/p`$ ratio with the beam energy. At mid-rapidity, the ratio asymptotes to unity with the increasing beam energy. We note the values of parameters shown in Table 1 are compatible within the statistical errors wit the results of Braun-Munzinger, Redlich, and Stachel , who obtain $`T_{\mathrm{chem}}=177\pm 7`$ MeV and $`\mu _B=29\pm 6`$ MeV.

## 4 Spectra

Let us turn now to the study of the transverse-momentum spectra of hadrons. They are calculated with the help of the Cooper-Frye  formula,

$$\frac{dN}{d^2p_{}dy}=p^\mu 𝑑\mathrm{\Sigma }_\mu f\left(pu\right),$$
(8)
where the distribution function $`f`$ includes the products of the resonance decays (more details on the construction of $`f`$ in the presence of the resonance decays have been given in Refs. ). The element of the hypersurface, $`d\mathrm{\Sigma }_\mu `$, is defined as

$$d\mathrm{\Sigma }_\mu =ϵ_{\mu \alpha \beta \gamma }\frac{x^\alpha }{\alpha _{}}\frac{x^\beta }{\alpha _{}}\frac{x^\gamma }{\varphi }d\alpha _{}d\alpha _{}d\varphi ,$$
(9)
where $`x^0=t`$, $`x^1=r_x`$, $`x^2=r_y`$, $`x^3=r_z`$, and $`ϵ_{\mu \alpha \beta \gamma }`$ is the Levi-Civita tensor. Introducing a convenient parameterization :

$`t`$ $`=`$ $`\tau \mathrm{cosh}\alpha _{}\mathrm{cosh}\alpha _{},r_z=\tau \mathrm{sinh}\alpha _{}\mathrm{cosh}\alpha _{},`$
$`r_x`$ $`=`$ $`\tau \mathrm{sinh}\alpha _{}\mathrm{cos}\varphi ,r_y=\tau \mathrm{sinh}\alpha _{}\mathrm{sin}\varphi ,`$ (10)
one finds

$$d\mathrm{\Sigma }^\mu (x)=u^\mu (x)\tau ^3\mathrm{sinh}(\alpha _{})\mathrm{cosh}(\alpha _{})d\alpha _{}d\alpha _{}d\varphi ,$$
(11)
hence the four-vectors $`d\mathrm{\Sigma }^\mu `$ and $`u^\mu `$ are parallel. This feature is special for our modeling of the freeze-out, see Eqs. (4) and (7), and allows us to represent the spectra in a compact form given by Eq. (8). With the use of the parameterization (10) we further rewrite Eq. (8) in the form

$$\frac{dN}{d^2p_{}dy}=\tau ^3\underset{\mathrm{}}{\overset{+\mathrm{}}{}}𝑑\alpha _{}\underset{0}{\overset{\rho _{\mathrm{max}}/\tau }{}}\mathrm{sinh}\alpha _{}d\left(\mathrm{sinh}\alpha _{}\right)\underset{0}{\overset{2\pi }{}}𝑑\xi puf\left(pu\right),$$
(12)
where

$$pu=m_{}\mathrm{cosh}\alpha _{}\mathrm{cosh}\alpha _{}p_{}\mathrm{cos}\xi \mathrm{sinh}\alpha _{}.$$
(13)
We note that $`\alpha _{}`$ is the rapidity of the fluid element, $`v_z=r_z/t=\mathrm{tanh}\alpha _{}`$, and $`\alpha _{}`$ describes the transverse size of the system, $`\rho =\tau \mathrm{sinh}\alpha _{}`$.

The maximum and average transverse-flow parameter is given in our model by the equations

$`\beta _{}^{\mathrm{max}}={\displaystyle \frac{\rho _{\mathrm{max}}}{\sqrt{\tau ^2+\rho _{\mathrm{max}}^2}}},\beta _{}={\displaystyle \frac{_0^{\rho _{\mathrm{max}}}\rho 𝑑\rho \frac{\rho }{\sqrt{\tau ^2+\rho ^2}}}{_0^{\rho _{\mathrm{max}}}\rho 𝑑\rho }}.`$ (14)
Since the ratios of hadron abundances depend very weakly on the centrality of the collision, we regard our thermodynamic parameters $`T`$ and $`\mu _B`$ as universal. On the other hand, the two geometric parameters may be different for each centrality class. From the obvious reasons we expect that $`\tau `$ and $`\rho _{\mathrm{max}}`$ are larger for more central events or, in other words, they are decreasing functions of the centrality, $`c`$, defined by the formula

$$c=\frac{\pi b^2}{\sigma _{\mathrm{tot}}},$$
(15)
where $`b`$ is the impact parameter and $`\sigma _{\mathrm{tot}}`$ is the total inelastic cross section for the colliding gold nuclei (for a discussion of the validity of this simple geometric relation see Ref. ). For various centrality bins we fit the geometric parameters with the help of the $`\chi ^2`$ method, i.e. we minimize the expression

$$\chi ^2(\tau ,\rho _{\mathrm{max}})=\underset{n=1}{\overset{n_{\mathrm{max}}}{}}\frac{\left[D_n^{\mathrm{exp}}D_n^{\mathrm{model}}(\tau ,\rho _{\mathrm{max}})\right]^2}{\sigma _n^2},$$
(16)
where $`D_n^{\mathrm{exp}}`$ is the n-th measured value of the transverse-momentum spectrum, $`D_n^{\mathrm{model}}(\tau ,\rho _{\mathrm{max}})`$ is the corresponding value calculated in our model, and $`\sigma _n`$ is the error. Observing the differences between the points measured by different experimental groups we conclude that the systematic errors are much larger than the statistical errors. Therefore, we assume that $`\sigma _n`$ is determined by the systematic uncertainties and take $`\sigma _n=D_n^{\mathrm{exp}}/10`$.444 We note that the $`\chi ^2`$ method was also used to determine the values of the two thermodynamic parameters in Sect. 3. Since the thermodynamic parameters are sensitive only to the ratios of hadron abundances our fitting of the four parameters $`T`$, $`\mu _B`$, $`\tau `$, and $`\rho _{\mathrm{max}}`$ can be done in two steps: first we determine $`T`$ and $`\mu _B`$, and having fixed these two values we determine later $`\tau `$ and $`\rho _{\mathrm{max}}`$. One could follow also a different path by fitting simultaneously the four parameters to the experimental spectra.

## 5 Most central events

In Fig 1 we show the result of our fit to the transverse-momentum spectra of $`\pi ^{},\overline{p}`$ and $`K^{}`$, for the most central events measured by the four RHIC experiments 555Most of the data shown in our figures were digitized from the available plots with preliminary experimental results. (exceptionally, the data of PHOBOS are for the mixtures of negative and positive hadrons). In this case, the centrality class of the BRAHMS and PHOBOS data is 10%, whereas the centrality class of the PHENIX and STAR data is 5%. Our model curves for $`\pi ^{},\overline{p}`$ and $`K^{}`$, shown in Fig. 1, are simultaneous fits to the four experimental spectra available for each hadron species. The fitted values of the geometric parameters are in this case $`\tau `$ = 7.58 fm and $`\rho _{\mathrm{max}}`$ = 7.27 fm (see Table 2 for a complete list of the values of the geometric parameters found for different experimental situations). We observe that these numbers are similar to those found in our previous study of the RHIC data collected at a lower energy of $`\sqrt{s_{NN}}`$ = 130 GeV, $`\tau `$ = 7.66 fm and $`\rho _{\mathrm{max}}`$ = 6.69 fm . This weak dependence of our geometric parameters on the energy of the colliding ions is reminiscent of the very weak energy dependence of the measured transverse HBT radii $`R_{\mathrm{side}}`$ and $`R_{\mathrm{out}}`$ .

In Fig. 2 we show our fit to the most central BRAHMS data only . In this case the fitted values of the geometric parameters are $`\tau =7.68`$ fm and $`\rho _{\mathrm{max}}=7.46`$ fm. These values are consistent with the previous values obtained from the “mixed” input. In Figs. 1 and 2 we can see that our model curves describe the data very well in the large range of $`p_{}`$ stretching up to about 2 GeV.

A closer inspection of Figs. 1 and 2 shows that our model curve tends to overestimate the spectra of antiprotons at small values of $`p_{}`$. To study the possible origin of this effect, in Fig. 3 we show in detail the antiproton spectra measured by different experimental groups. One can see that the experimental spectra differ substantially from each other and, in practice, one cannot fit all four groups of the experimental points simultaneously. The most likely reason for such differences are difficulties connected with the correct estimate of the role of the weak decays, or with fixing the normalization. Our model calculations shown in Fig. 3 indicate that the effect of the weak decays on the antiproton spectrum is indeed important. The solid line represents our result obtained in the case where the antiprotons are fed by the weak decays. This curve may be treated as the upper limit for the spectrum. On the other hand, the dashed line represents our result with no feeding included, and may be treated as the lower limit. We can see that most of the experimental points, except for those of PHOBOS, lie in the band enclosed with our two curves.

## 6 Centrality dependence of the geometric parameters

Having at our disposal the data collected in different centrality windows we may determine the dependence of the expansion parameters on $`c`$, or, equivalently, on the impact parameter $`b`$. In Fig. 4 we show the PHENIX data together with our model curves. Both the data and the model spectra are plotted for eight different centrality classes. The most central collision correspond to the centrality class $`c=`$ 0-5%, whereas the most peripheral collisions included in Fig. 4 are for the centrality class $`c=`$ 80-91%. The values of the fitted expansion parameters for all eight cases are given in Table 2.

The geometric parameters obtained from the fit to the STAR data turn out to be noticeably larger compared to the parameters obtained from the analysis of the PHENIX data. In order to see the origin of this effect, in Fig. 5 we superimpose the spectra from STAR and PHENIX for four centrality classes common for the two experiments. We can see that the spectra measured by STAR (open symbols) are lying slightly above the spectra measured by PHENIX (filled symbols). As a consequence of this behavior, our normalization-controlling parameter, $`\tau `$, turns out to be larger for STAR. We have performed a simultaneous fit to the PHENIX and STAR data. The model curves for this case are also shown in Fig. 5 (solid, dashed, and long-dashed curves, respectively for pions, kaons, and antiprotons). The fitted values of the expansion parameters $`\tau `$ and $`\rho _{\mathrm{max}}`$ are listed in the lowest entry of Table 2.

An interesting behavior can be found when the values of $`\tau `$ and $`\rho _{\mathrm{max}}`$ are plotted versus $`c`$ (see Figs. 6 and 7). The parameter $`\rho _{\mathrm{max}}`$ exhibits almost ideal linear dependence on $`c`$ in the full range $`0<c<1`$. Similarly, the parameter $`\tau `$ decreases linearly with growing $`c`$ in the range $`0<c<0.5`$. For larger values of $`c`$ the drop of $`\tau `$ is slightly weaker, although the data of PHENIX show a rather strong tendency to line up. The approximate linear scaling of both $`\tau `$ and $`\rho _{\mathrm{max}}`$ with $`c`$ has an important physical significance. The hadron multiplicities in our model are functions of $`\tau ^3`$ and $`\rho _{\mathrm{max}}/\tau `$, see Eq. (12), hence the linear dependence of $`\tau `$ and $`\rho _{\mathrm{max}}`$ on $`c`$ means that the ratio $`\rho _{\mathrm{max}}/\tau `$ is independent of $`c`$ and we are left with the dependence on $`\tau ^3`$ only. Consequently, using Eq. (15) we find that the total multiplicity in the central region (mid-rapidity) obtained in our model should, to a good approximation, scales in the following way

$$N_{\mathrm{model}}(b)=N_{\mathrm{model}}(b=0)\left(1c(b)\right)^3=N_{\mathrm{model}}(b=0)\left(1\frac{\pi b^2}{\sigma _{\mathrm{tot}}}\right)^3.$$
(17)
This type of the behavior may be compared to the $`b`$-dependence of hadron production following from the wounded nucleon model . For the symmetric case, i.e. when two identical nuclei collide, the number of the wounded nucleons is given by the formula

$$w\left(b\right)=2Ad^2sT_A\left(𝐛𝐬\right)\left(1\left[1\sigma _{\mathrm{in}}T_A\left(𝐬\right)\right]^A\right)$$
(18)
where $`A`$ is the mass number, $`\sigma _{\mathrm{in}}`$ is the inelastic nucleon-nucleon cross section, and $`T_A`$ is the nucleon-nucleus thickness function. Similarly, the total inelastic cross section is given by the formula

$$\sigma _{\mathrm{tot}}=d^2b\left(1\left[1\sigma _{\mathrm{in}}T_{AA}\left(𝐛\right)\right]^{A^2}\right),$$
(19)
where $`T_{AA}\left(𝐛\right)`$ is the nucleus-nucleus thickness function. Using the values $`A=197`$ and $`\sigma _{\mathrm{in}}`$=40 mb, and assuming the Woods-Saxon distribution for the nuclear densities (with the standard choice of the parameters) we find

$$w\left(b\right)w(0)(1c(b))^3.$$
(20)
Since the wounded-nucleon model assumes the scaling of hadron production with the number of the wounded nucleons, the total multiplicity of hadrons considered in this model scales also as $`(1c(b))^3`$, hence it is compatible with the scaling obtained in our approach, cf. Figs. 6 and 7. The quality of the simple approximation (20) is shown in Fig. 8. This result means that our model, with the geometric parameters fit to reproduce the data at various centralities, conforms with high accuracy to the wounded-nucleon scaling for not-too-large values of $`c`$ (say, $`c<0.5`$).

We end this Section with a technical remark. In the analysis of this Section we have ignored any possible effects of the azimuthal deformation. Clearly, the emitting source formed in a non-central collision will not be azimuthally-symmetric, which, by the way, is seen in the HBT measurements . However, it turns out that the effects of this deformation are tiny for the $`p_T`$-spectra averaged over the azimuthal angle , thus can be dropped in the analysis of the present paper.

## 7 Spectrum of the $`\mathrm{\Omega }`$ baryons

Finally, we wish to present the $`p_{}`$ spectra for the $`\mathrm{\Omega }`$ and $`\overline{\mathrm{\Omega }}`$. Similarly to the case of the lower energy of $`\sqrt{s_{NN}}=130`$ GeV , our model is in full agreement with the preliminary data of Ref. , which can be seen in Fig. 9. In the model calculation we have used the thermal parameters from Table 1 and the geometric parameters $`\tau =8.3`$ fm and $`\rho _{\mathrm{max}}=7.1`$ fm, which roughly corresponds to centrality $`010`$% in the STAR experiment (cf. Table 2).

## 8 Conclusions

It is clear from the results presented in this paper that the single freeze-out model works well for the particle ratios and the transverse-momentum spectra for the full RHIC energy of $`\sqrt{s_{NN}}=200`$ GeV. This confirms our results published previously for $`\sqrt{s_{NN}}=130`$ GeV . The model provides an economic framework to parameterize the data, and supports the use of the thermal approach to heavy-ion physics. In addition, we have shown that the model works equally well for the central, as well as for the non-central collisions. The transverse-size parameter, $`\rho _{\mathrm{max}}`$, drops linearly practically for all values of $`c`$, while the ratio of $`\rho _{\mathrm{max}}`$ to $`\tau `$ remains almost constant for not-too-large values of $`c`$. The systematics of the dependence of the geometric/flow parameters on the centrality is most interesting. For moderate values of the parameter $`c`$, say $`c<0.5`$, the analysis of the non-central collisions displays scaling in accordance to the wounded-nucleon model.

## Acknowledgements

We are grateful to Fuqiang Wang and Patricia Fachini for many helpful remarks concerning the data.


text:
# Theoretical analysis of continuously driven dissipative solid-state qubits

\[

## Abstract

We study a realistic model for driven qubits using the numerical solution of the Bloch-Redfield equation as well as analytical approximations using a high-frequency scheme. Unlike in idealized rotating-wave models suitable for NMR or quantum optics, we study a driving term which neither is orthogonal to the static term nor leaves the adiabatic energy value constant. We investigate the underlying dynamics and analyze the spectroscopy peaks obtained in recent experiments. We show, that unlike in the rotating-wave case, this system exhibits nonlinear driving effects. We study the width of spectroscopy peaks and show, how a full analysis of the parameters of the system can be performed by comparing the first and second resonance. We outline the limitations of the NMR linewidth formula at low temperature and show, that spectrocopic peaks experience a strong shift which goes much beyond the Bloch-Siegert shift of the Eigenfrequency.

\]

Coherent manipulation of quantum states is a well established technique in atomic and molecular physics. In these fields, one works with “clean” generic quantum systems which can be very well decoupled from their environments. Moreover, it is possible to apply external fields in a way such that strong symmetry relations between the static and the time-dependent part of the Hamiltonian apply and the resulting dynamics is very simple and can be treated analytically. In solid-state systems, the situation is different. Not only do they contain a macroscopic number of degrees of freedom which form a heat bath decohering the quantum states to be controlled, but also is the choice of controllable parameters much more restricted. A quantum-mechanical two state system (TSS) realized in a mesoscopic circuit can be identified with a (pseudo)spin, however, in that case the different components of the spin may correspond to physically distinct observables such as e.g. magnetic flux and electric charge . This naturally limits the possibilities of controlling arbitrary parameters of the pseudospin. Hence, in order to describe the direct control of quantum states in mesoscopic devices, concepts from NMR or quantum optics cannot be directly applied but have to be carefully adapted. In particular, as decoherence is usually rather strong in condensed matter systems, one can attempt to drive the system rather strongly in order to have the operation time for a quantum gate, usually set by the Rabi frequency, as short as possible.

We concentrate on the case of a persistent current quantum bit driven through the magnetic flux through the loop and damped predominantly by flux noise  with Gaussian statistics. This setup is accurately described by the driven spin-boson model

$$H=\frac{ϵ(t)}{2}\widehat{\sigma }_z\frac{\mathrm{\Delta }}{2}\widehat{\sigma }_x+\widehat{\sigma }_z\underset{i}{}c_ix_i+\underset{i}{}\left(\frac{\widehat{p}_i^2}{2m_i}+\frac{1}{2}m_i\omega _i^2x_i^2\right)$$
(1)
where $`ϵ(t)=ϵ_0+s\mathrm{cos}\mathrm{\Omega }t`$ and the oscillator bath is assumed to be ohmic with a spectral density $`J(\omega )=\frac{\pi }{2}_i\frac{c_i^2}{m_i\omega _i}\delta (\omega \omega _i)=\alpha \omega e^{\omega /\omega _c}`$. The connection of $`J(\omega )`$ to the setup parameters is detailed in . The static energy splitting of the peudospin is $`\nu =\sqrt{ϵ_0^2+\mathrm{\Delta }^2}`$. This model is also applicable to other Josephson qubits and other realizations . In particular, the strong driving regime we are going to elaborate on has recently been realized in several setups . We study the effective dynamics of the pseudospin having traced out the bath in the limit of weak damping, $`\alpha 1`$, which is appropriate for quantum computation. This is done using the Bloch-Redfield equation. The resulting equation is of Markovian form in the sense that it only contains the density matrix at a single time, however, it is derived in such a way, that the free coherent evolution during the interaction with the bath is fully taken into account such that the resulting equation is numerically equivalent to a fully non-Markovian path-integral scheme  and only memory terms beyond the Born approximation are dropped. The explicit form of the equations for this situation as well as the formulas for the rates correspond to those given in . We compare our numerical results to analytical formulae derived in the framework of a high-frequency approximation which involves averaging over the driving field and has nonetheless shown to give a good estimate for the system dynamics even close to resonances.

Initial experiments on quantum bits such as Ref.  do not monitor the real-time dynamics of the system as in Ref. , because the read-out is much slower than the decoherence, i.e. the dephasing time $`\tau _\varphi `$ is too short. In order to optimize the experimental setup, it is important to measure both $`\tau _\varphi `$ and the relaxation time $`\tau _R`$, even and in particular if they are insufficient. In the standard NMR-case, this is done by studying the width of the resonance . We will detail, a somewhat modified analysis can be performed for solid-state qubits and what are its limitations. We discuss both situations. Our results thus help to analyze the decoherence as observed in Refs. , and outline the possibilities and limitations of driving the system in the nonlinear regime.

We have numerically solved the driven Bloch-Redfield equation. The real-time dynamics is illustrated in fig. 1.

The dynamics shows distinct features on different time scales. As expected, there are clear Rabi oscillations on the scale of the effective driving strength (see below). In quantum computing applications, these would be used for the implementation of a Hadamard gate. On top of this, there are fast components: The dominating one oscillates with the driving frequency, which originates in the fact that the driving is not perpendicular to the static field. A weaker one, which oscillates at twice the driving frequency, comes from the counter-rotating term perpendicular to the static field. These oscillations can lead to errors of the Hadamard gate. On a longer time scale, the Rabi oscillations decay. The time scales will be discussed later on. In general, if one is not exactly on resonance, these oscillations are combined with nonoscillatory decay, see fig. 1 a) and b). At very long times, the system assumes a quasistationary value $`P_{\mathrm{}}`$.

Corresponding to the situation of a spectroscopy experiment, we now turn to the analysis of the quasi-stationary state which is established after a long time $`t\tau _\varphi ,\tau _R,\omega _R^1`$. We compare our full numerical solutions with analytical expressions we have obtained from the high-frequency approximation of Refs. . As a result of this approach, the TSS is mapped onto a coupled ensemble of TSSs corresponding to the original system emitting or absorbing $`n`$ photons from the driving field during the tunneling. The energy bias of these individual systems is $`ϵ_n=ϵ_0n\mathrm{\Omega }`$ and the tunnel matrix element

$$\mathrm{\Delta }_n=\mathrm{\Delta }J_n(s/\mathrm{\Omega })$$
(2)
where the $`J_n`$ are Bessel functions. At low driving fields, we can approximate $`\mathrm{\Delta }_n=\frac{\mathrm{\Delta }}{n!}(s/2\mathrm{\Omega })^n`$ as we would expect from the expansion of a perturbation series in the driving strength. The $`\mathrm{\Delta }_n`$ can hence be viewed as $`n`$-photon Rabi frequencies. This implies, that the usual single-photon frequency gets replaced by $`\mathrm{\Delta }_1s\mathrm{\Delta }/\nu `$, which can be interpreted as only the projection of the driving field onto the direction in pseudospin space orthogonal to the static Hamiltonian. In order to obtain the solid curves in fig. 2 the secular equations for the Eigenfrequencies have been solved, taking into account an appropriate number of terms . The dynamical two state systems are characterized by individual dynamical dephasing rates $`\mathrm{\Gamma }_{\varphi ,n}`$ and a common relaxation rate $`\mathrm{\Gamma }_r`$ . On the n-th resonance, $`\mathrm{\Gamma }_{\varphi ,n}`$ can be very low, much lower than off-resonance, as can be seen in figure 1, and largely exceed the intrinsic dephasing time. This has been observed in Refs. .

Figure 2 shows numerical and analytical results for $`P_{\mathrm{}}=lim_t\mathrm{}\sigma _z`$ at a fixed frequency $`\mathrm{\Omega }/2\pi =6.6`$GHz as a function of the energy bias $`ϵ_0`$. This corresponds to a realistic experimental situation . In fig. 2 a), taken at weak driving field, only the regular resonance corresponding to the transition between the two Eigenstates driven by absorbing a single photon can be seen. At somewhat stronger driving, fig. 2 b), this peak grows wider and a second resonance appears, corresponding to the simultaneous absorption of two photons. At higher fields, fig. 2 c), these peaks grow and start to dominate over the background. They also turn asymmetric. This trend culminates in the situation shown in fig. 2 d). In that case, $`P_{\mathrm{}}`$ does not grow to positive values at small positive $`ϵ_0`$, but it gets negative and then directly approaches the first resonance. The reason for this behavior can be identified within the high-frequency approximation: The lowest order-tunnel frequency $`\mathrm{\Delta }_0=\mathrm{\Delta }J_0(s/\mathrm{\Omega })`$ vanishes at this particular driving strength. Indeed, comparing figs. 2 a) and d) one can see, that the step which is at $`ϵ_0=0`$ in case a) is shifted to $`ϵ_0\mathrm{\Omega }`$ in case d). This phenomenon, the coherent desctruction of tunneling  relies on destructive interference of the dressed state  formed by the TSS and a cloud of photons from the driving field. This interpretation is supported by the dynamics of $`\sigma _z(t)`$. As seen in fig. 3, which shows the dynamics at the degeneracy point for different driving strengths, the zero-photon tunneling is slowed down and brought to a standstill. If that strong driving can be applied to solid-state qubits, it would provide an alternative for controlling $`\mathrm{\Delta }_0`$ by a cw microwave field instead of an additional magnetic flux as proposed in Ref. .

At very weak driving, the peak position corresponds to the qubit eigenfrequency, $`\mathrm{\Omega }=\nu `$. This is not reliably predicted by the high-frequency approximation. At stronger driving, the peak gets shifted. Closer inspection as in figure 4 shows, that this shift goes much beyond the usual Bloch-Siegert shift  of the dynamical Eigenfrequency, in fact, one can show that the position of the peak in steady state and the Eigenfrequency do not conincide. The former is given by balancing of rates and it can be shown that in lowest order gets shifted by $`\delta ϵ_{\mathrm{peak}}s^2/8\mathrm{\Omega }^2`$ whereas the Bloch-Siegert shift for our case is $`\delta ϵ_{\mathrm{BS}}\mathrm{\Delta }^2s^2/(16\mathrm{\Omega }^3)`$. As a more general conclusion, already at modest not-too-weak driving, the resonance positions do not necessarily reflect the Eigenfrequencies of the system.

In fig. 5, the height of the two lowest order peaks is shown. It can be seen, that, from the low-driving side, they saturate as soon as their effective Rabi frequency $`\mathrm{\Delta }_n`$ exceeds $`1/\sqrt{T_rT_\varphi }`$. At very high driving, the peaks show an inversion of population.

For the optimization of qubit setups on the way to coherent dynamics, it is important to characterize its coherence properties from the spectroscopic data. In NMR, this is done from the linewidth given by

$$\delta \mathrm{\Omega }=2\sqrt{\tau _\varphi ^2+\omega _R^2\tau _R/\tau _\varphi }$$
(3)
where $`\omega _R`$ is the Rabi frequency at resonance, which coincides with the strength of the driving field. A generalization of this formula to our case has to take into account low temperatures and the different driving situation. Moreover, $`\omega _R`$ is usually not directly known to sufficient precision, because the driving strength depends on the attenuation of the applied fields on their way to the sample and the efficiency of the coupling .

Our analysis suggests the generalization of eq. 3 is given by

$$\delta \mathrm{\Omega }_n=2\sqrt{\tau _\varphi ^2+\mathrm{\Delta }_n^2\tau _R/\tau _\varphi }n=\pm 1,\pm 2\mathrm{}$$
(4)
where $`\delta \mathrm{\Omega }_n`$ is the width (in frequency) of the $`n`$-photon resonance and $`\mathrm{\Delta }_n`$ is the effective Rabi frequency defined above. At low powers $`s<\mathrm{\Omega }`$, they are given by the rates from the undriven ohmic case

$$\tau _\mathrm{R}^1=\alpha \frac{\mathrm{\Delta }^2}{2\nu }\mathrm{coth}\left(\frac{\mathrm{}\nu }{2k_BT}\right)\tau _\varphi ^1=(2\tau _R)^1+2\pi \alpha \frac{k_BT}{\mathrm{}}\frac{ϵ_0^2}{\nu ^2}.$$
(5)
This result is confirmed by our numerical simulations fig. 6. We can essentially identify three regimes: A saturation broadening regime at low powers, where $`\delta \mathrm{\Omega }_n2\tau _\varphi ^1`$, a saturated regime, $`\delta \mathrm{\Omega }_n2\mathrm{\Delta }_n\sqrt{\tau _R/\tau _\varphi }`$ and a nonlinear regime, where the numerical curve deviates from eq. 4 due to the fact, that the high Rabi frequency shifts the relevant energy scales and modifies the timescales given in eq. 5. Note, that in this regime, the general curve of $`P_{\mathrm{}}`$ is greatly deformed (see fig. 2) and the width of a peak becomes ambiguous.

This result allows to measure essentially all interesting parameters of the system experimentally. By extrapolating the level separation at the degeneracy point (as it was done in ), one obtains $`\mathrm{\Delta }`$. By tracking the resonance positions at weak driving, one can evaluate $`ϵ_0`$ as a function of the external control parameter (in  this would be the magnetic flux). By driving in the saturated regime, the widths of the first and second peak become, according to eq. 4 $`\delta \mathrm{\Omega }_{1/2}=2\mathrm{\Delta }_{1/2}\sqrt{\tau _R/\tau _\varphi }`$, hence by taking their ratio we find the effective driving strength from $`\mathrm{\Delta }_2/\mathrm{\Delta }_1=J_1(s/\mathrm{\Omega })/J_0(s/\mathrm{\Omega })s/2`$ and by tracking the slope of the first resonance we find the ratio $`\tau _r/\tau _\varphi `$. Finally, examining the saturation broadening regime of the first resonance gives the absolute value of $`\tau _\varphi `$.

In conclusion, we have numerically and analytically analyzed the spin-boson system, which e.g. represents a SQUID-qubit, in the weak damping regime, driven by continuous fields. As compared to the more familiar situation in NMR, this system is both different in the character of the driving and the low temperature governing the dissipation. We have shown, that the key features of this system, Rabi-oscillations, and saturation of the linewidth, persist qualitatively as has been experimentally confirmed . They are however altered on a quantitative level, such as an unanticipatedly strong shift of the position of the resonance peak, and also supplemented by new phenomena such as higher-harmonics generation, oscillations of $`\sigma _z`$ on the scale of the driving field, and coherent desctruction of tunneling. We have finally outlined a scheme how to determine all relevant parameters (tunnel splitting, energy dispersion, driving strength, dephasing and relaxation time) of a quantum bit solely through spectroscopy.

We thank M. Grifoni, C.H. van der Wal, A.C.J. ter Haar, C.J.P.M. Harmans, J. von Delft and I. Goychuk for discussions. Work supported by the EU through TMR “Supnan” and IST “Squbit”. FKW acknowledges support by the ARO under contract-Nr. P-43385-PH-QC.


text:
# The Massive Triple Star System HD 16429 A

## 1 Introduction

O-type stars are most commonly found within young open clusters and OB associations in the spiral arms of the Galaxy (Gies, 1987), and the binary frequency among these cluster and association members is significantly higher (between $`7075\%`$) than among field O stars ($`1922\%`$) and runaways (about $`8\%`$) (Gies, 1987; Garmany, 1994; Mason et al., 1998). Multiple O star systems in clusters are particularly influential in the dynamics of the cluster because stars will be ejected both through supernova explosions in close binaries and through close gravitational encounters of stars with binaries. Triple O star systems containing a close pair and a more widely separated component are particularly important sites of gravitational interactions (Mason et al., 1998).

HD 16429 A (HIP 12495, V482 Cas, SAO 12383, BD+$`60^{}541`$) is an O9.5 II ((n)) star (Walborn, 1976) and a member of the Cas OB6 association (Garmany & Stencel, 1992). It has two optical companions, component B at a separation of 6.66 arcsec, and C at a separation of 53.2 arcsec (Mason et al., 1998). HD 16429 B is an F4 V star (Meisel, 1968), which indicates that it is not physically associated with A. Component C is more likely associated based upon its magnitude and color (see §4). The A component was itself found to be a speckle binary with a separation of 0.30 arcsec (Mason et al., 1998). Gies & Bolton (1986) found that one component of A may also be a spectroscopic binary, making it a triple system, but they were unable to identify the period of this close pair. All of the components of the heirarchical HD 16429 system are shown in Figure 1. Other observations have shown HD 16429 A to be a radio emitter (Martí et al., 1998) as well as a weak X-ray source (Goldoni & Mereghetti, 1995). The radio and X-ray emission probably originates in the stellar winds (Martí et al., 1998).

Here I present evidence to confirm that HD 16429 A is a triple star system with a spectroscopic binary component (HD 16429 Ab1) and a radial velocity stationary companion (HD 16429 Aa). The Aa and Ab1 components were isolated using a Doppler tomography algorithm, but the Ab2 component remains hidden. Approximate spectral types and relative flux contributions are presented for the two visible components. Finally, I discuss the stars in the immediate vicinity of HD 16429 A and their relationship to the Cas OB6 association.

## 2 Observations and Blended Line Radial Velocities

Thirty-five optical spectra of HD 16429 A were obtained with the Kitt Peak National Observatory 0.9-m Coudé Feed Telescope during runs in 2000 October and 2000 December. The spectra have a resolving power $`R=\lambda /\delta \lambda =9500`$ using the long collimator, grating B (in second order with order sorting filter OG550), camera 5, and the F3KB CCD, a Ford Aerospace $`3072\times 1024`$ device. These arrangements produced a spectral coverage of $`64407105`$ Å. Usually two exposures of 15 minutes duration were taken each night, separated by $`23`$ hours. The spectrum of the rapidly rotating A-type star, $`\zeta `$ Aql, was observed each night for removal of atmospheric water vapor and O2 bands. Also, the O-type stars HD 34078 and HD 47839 were observed as spectral comparisons. Each set of observations was accompanied by numerous bias, flat field, and Th Ar comparison lamp calibration frames. The dates of observation are given in Table 1.

The spectra were extracted and calibrated using standard routines in IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which is operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation.. All the spectra were rectified to a unit continuum by fitting line-free regions. The removal of atmospheric lines was done by creating a library of $`\zeta `$ Aql spectra from each run, removing the broad stellar features from these, and then dividing each target spectrum by the modified atmospheric spectrum that most closely matched the target spectrum in a selected region dominated by atmospheric absorptions. Some glitches appeared as artifacts of the removal of the atmospheric telluric lines, and these were excised by linear interpolation. The spectra from each run were then transformed to a common heliocentric wavelength grid.

Radial velocities were measured by fitting a parabola to the line cores of H$`\alpha `$, He I $`\lambda 6678`$, and He I $`\lambda 7065`$, but several factors made it difficult to obtain consistent measurements. The H$`\alpha `$ line shows variable P Cygni emission that pushes the line’s center blueward. Due to this certain contamination, radial velocities from the H$`\alpha `$ line were not used. In addition, measurements of the He I $`\lambda 7065`$ line have a slightly lower mean radial velocity than the mean for the He I $`\lambda 6678`$ line. This could be due to different relative line strengths of the He I lines in the blended stars, contamination from He II $`\lambda 6683`$ absorption in one component, or a radial velocity gradient due to outflowing winds (see §4). However, the difference in their mean velocities is small and only corresponds to a difference of 1.4 pixels on the chip, so the mean velocity of the two He I lines was used in this analysis.

To account for small deviations in radial velocity due to the wavelength calibration, the interstellar lines present in the spectra were inspected. A large number of interstellar features were identified from the lists of Herbig (1975), Herbig & Leka (1991), Morton (1991), and especially Galazutdinov et al. (2000). An interstellar spectrum was created by extracting these lines from the mean HD 16429 spectrum, and then it was cross correlated with each individual spectrum to measure any systematic radial velocity shift (generally $`<7`$ km s-1). The velocity shift was subtracted from the measured radial velocities, and the resulting velocities are shown in the third column of Table 1. In the next section, I will show that these radial velocities actually represent the results of line blending between a single-lined spectroscopic binary and a stationary component.

## 3 Orbital Elements from a Two Component Model

Visual inspection of the spectra revealed radial velocity variations indicative of an orbital period of about 3 days, but a more precise period search was done using a version of the discrete Fourier transform and CLEAN deconvolution algorithm of Roberts, Lehár, & Dreher (1987) (written in IDL222IDL is a registered trademark of Research Systems, Inc. by A. W. Fullerton). The most significant signal found in the CLEANed power spectrum occurs at the period $`P=3.051`$ d. Relics of this signal also occur at integer multiples of this frequency, but the true period cannot be shorter because there are no significant changes in radial velocity over the $`23`$ hour intervals between the observations from a given night.

This period was used with the non-linear, least-squares fitting program of Morbey & Brosterhus (1974) to solve for the orbital elements. The best fit was obtained by allowing the period to vary slightly, to an improved value of $`P=3.0539`$ d, and assigning only half weight to one point with very large scatter near 0.95 phase. The resulting preliminary orbital elements are listed in the second column of Table 2, and the radial velocity curve is illustrated in Figure 2 with the dotted line. Note that this curve has been shifted in phase so that the epoch, $`T`$, corresponds to the final time of periastron, discussed below. The observed minus calculated velocity errors ($`OC`$) based on this preliminary solution are given in Table 1, column 4. The root-mean-square of the velocity residuals is comparable to the scatter in the best measured lines.

There are a number of other radial velocity measurements for HD 16429 A in the literature (10 by Abt, Levy, & Gandet (1972), 15 by Gies & Bolton (1986), and 1 by Conti, Leep, & Lorre (1977)), so these measurements were combined with those in Table 1 to refine the orbital period. Since each investigation used different spectral lines and wavelength regions to measure velocities, the results may have systematic differences. However, the agreement between data sets was reasonable, and including the additional measurements significantly improved the orbital period, $`P=3.05442\pm 0.00005`$ d.

The speckle results (Mason et al., 1998) indicate that at least two bright components were recorded in these spectra, so I created grayscale images of the spectral lines as a function of orbital phase to search for evidence of multiple components. Figure 3 shows the plot of the He I $`\lambda 6678`$ line. It clearly reveals signs of a single-lined spectroscopic binary, Ab1, superimposed upon a stationary component, Aa, with similar systemic velocity, consistent with expectations from the speckle observations. Because the stationary component’s lines are blended with the moving component, the preliminary radial velocity semiamplitude of the moving component is probably greatly underestimated, and perhaps the preliminary eccentricity and longitude of periastron are misrepresented as well.

To determine the true orbital elements of the Ab1 component, the Aa component was isolated using a two-component tomographic separation and then subtracted from the observations. Tomographic separation relies upon the radial velocity estimates of the stars and their flux ratio to recreate the spectra of the individual stars (Bagnuolo & Gies, 1991). The algorithm uses the observed composite spectra, approximate flux ratios of the stars, the center of mass velocity, $`V_0`$, and the velocity semi-amplitude, $`K`$, of component Ab1 to calculate the solution using an iterative least squares technique. For the Aa component, the velocity amplitude was set to 0. A three-component separation was not attempted because of the extreme faintness of the Ab2 spectroscopic component. Even the two-component separation was difficult because the true velocity amplitude, $`K`$, of Ab1 is unknown, and visual inspection of Figure 2 reveals that it could be as high as 200 km s-1. The tomographic separation is not very sensitive to the choice of the relative flux contribution of the components. Initially, an equal flux contribution was assumed, and the value was later refined according to the methods discussed in the following section.

I computed tomographic separations for a grid of assumed semi-amplitudes for $`30<K<200`$ km s-1. In order to identify the best fit value of $`K`$, I used the separated components to recreate a set of synthetic spectra for each observation which was then compared to the original spectra. The He I $`\lambda 6678`$ and He I $`\lambda 7065`$ line regions were used separately to compute the $`\chi ^2`$ error between the observed and recreated spectra. The errors associated with the He I $`\lambda 7065`$ line were slightly larger because that line is located in a region with many atmospheric telluric lines, so it was correspondingly noisier after the telluric removal process. The values of $`K`$ corresponding to the minimum error for each line were then averaged to obtain the optimal value of $`K=124\pm 25`$ km s-1.

Once the Aa component was isolated, I subtracted its profile from the observed spectra to attain profiles of the Ab1 component alone, and then I measured the radial velocities of Ab1 in the same way as before. These unblended radial velocities are listed in column 5 of Table 1. The new orbital elements of the spectroscopic binary were calculated using the same technique discussed above and are included in the last column of Table 2. Table 1, column 2 lists the orbital phase from periastron based on this final solution. Although the optimal $`K`$ is $`124\pm 25`$ km s-1, the orbital fitting program found a slightly higher value of $`K=136\pm 6`$ km s-1, but these estimates are consistent within errors. Figure 2 includes the final radial velocities for the moving component as the solid points and the theoretical radial velocity curve as the solid line, and the new grayscale plot of the isolated moving component is shown in Figure 4. The signal to noise ratio is too low to positively identify the anti-phase motion of the spectroscopic companion Ab2 in this plot.

## 4 Discussion

The flux ratio of the two visible components of HD 16429 A was determined using their observed speckle magnitudes, $`\delta m_V=0.5\pm 0.2`$ (Mason & Hartkopf 2001, private communication). The stationary star contributes a stronger signal in the spectra, and it is also the brighter component with 61 $`\pm 5\%`$ of the overall luminosity. This star is HD 16429 Aa described by Mason et al. (1998), whereas the visible moving component is HD 16429 Ab1. The flux ratio can also be estimated using the equivalent width of each component’s He I $`\lambda 6678`$ line. Conti (1974) found that in late O stars of all luminosity classes, the equivalent width of this line is roughly constant (but with a large scatter), and I adopted the method of Petrie (1939) to determine the flux ratio of the stars from their observed equivalent widths. Assuming an intensity ratio $`r=F_2/F_1`$ for the tomographic reconstruction, the resulting equivalent widths of the reconstructed moving and stationary component spectra are $`A_1`$ and $`A_2`$, respectively. The relative continuum intensities are then $`1/(1+r)`$ and $`r/(1+r)`$, respectively, and the absolute line strengths are $`s_1=A_1/(1+r)`$ and $`s_2=rA_2/(1+r)`$. Because the true strengths of the lines are the same in each star, the true intensity ratio is $`r_1=rA_2/A_1`$. By integrating over the He I $`\lambda 6678`$ line in each component and comparing their equivalent widths, the contribution of the stationary component was found to be about $`70\pm 7\%`$ of the system’s total luminosity. This value agrees reasonably well with that obtained from speckle observations, but is less accurate because of the intrinsic scatter in the He I $`\lambda 6678`$ widths, so I adopted the flux ratio from the speckle observations.

Figure 5 presents a plot of the reconstructed spectra, using the final $`K=136`$ km s-1 and normalized for their flux contribution. The figure also shows the Aa and Ab1 components compared to two stars of similar spectral type, HD 34078 and HD 47839. The comparison spectra have been smoothed using a Gaussian profile so that their FWHM and line depths more closely match those of the separated components. This smoothing does not significantly affect the equivalent widths of the lines. The spectra of the Aa and Ab1 components have also been smoothed slightly to reduce noise.

The tomographic separation of the spectra revealed that both visible components of HD 16429 A have strong H$`\alpha `$ and He I lines. The smoothed spectra of HD 34078 and HD 47839 have He I lines of similar strength, although stronger H$`\alpha `$ profiles. H$`\alpha `$ emission is common in luminous O stars (luminosity classes I, II, and III), so the difference in line strengths is likely due to emission in HD 16429 A. (HD 34078 and HD 47839 are both main sequence O stars with no emission.) The emission may not have been distributed properly between the two component stars during the tomographic separation because it varies on a non-orbital timescale. The original blended spectra showed striking P Cygni emission in the H$`\alpha `$ line, which is associated almost entirely with the Aa component in the reconstructions. This widely-separated member of the triple system is probably the more evolved of the group since it has such strong emission from stellar winds. Since it is the more luminous, the stationary component corresponds to the O9.5 II type that was assigned to the system by Walborn (1976). The ((n)) suffix used by Walborn to indicate slightly broadened lines probably reflects the line blending with the Ab1 component. The He I line strengths are similar to those in HD 34078, an O9.5 V star. The absence of significant He II provides further confirmation of the spectral type.

The separation also revealed weak He II $`\lambda \lambda 6527`$, 6683, and 6891 absorption in the Ab1 component. He II lines only appear in O spectra, disappearing entirely in the earliest B spectra. Their presence indicates that this star is hotter than the brighter Aa star. Comparing the strengths of these lines with those in HD 47839 reveals that the moving component is a slightly later, O8 spectral type, although the luminosity class cannot be specified by the comparison. However, the difference in speckle magnitudes is consistent with an O8 III-V classification for HD 16429 Ab1 (Howarth & Prinja, 1989).

The second spectroscopic companion, Ab2, remains hidden even after the tomographic separation (Fig. 4), but the statistical method of Mazeh & Goldberg (1992) suggests that it is a massive, early B-type star. For single-lined spectroscopic binaries consisting of a primary of mass $`M_1`$ and a secondary of mass $`M_2`$, with a mass ratio $`q=M_2/M_1`$, the mass function is given by $`f(M_2)=M_1`$ sin $`{}_{}{}^{3}i`$ $`q^3/(1+q)^2`$. It cannot be solved uniquely for both $`q`$ and $`i`$, but Mazeh & Goldberg recommend an iterative solution based on a random distribution of both parameters. Given the known mass function and the primary mass, the mass ratio can be calculated for each possible inclination and then averaged over the range of inclinations. Using a mass of $`2833M_{\mathrm{}}`$ for the O8 III-V primary (Howarth & Prinja, 1989) and $`f(M_2)=0.76\pm 0.11`$ from the final orbital solution, the probable mass ratio is 0.6. This is inconsistent with a low mass compact companion, but suggests instead a $`1720M_{\mathrm{}}`$ early B-type star (Howarth & Prinja, 1989). The Ab2 component remains hidden because its flux is overwhelmed by the two brighter O stars. If these masses are adopted, then the orbital inclination is approximately $`i=42^{}`$, and no eclipses are expected (or observed by Hipparcos; Perryman et al. 1997).

There are several stars in the vicinity of HD 16429 A, including the Be X-ray binary LS I +$`61^{}303`$, that have similar radial velocities, proper motions, apparent $`V`$ magnitudes, and $`BV`$ colors (summarized in Table 3). Parabolic fits of the He I $`\lambda 6678`$ and He I $`\lambda 7065`$ lines in the tomographically isolated spectrum of the Aa component provided its average radial velocity of $`26\pm 5`$ km s-1. This is somewhat different than the systemic velocity of the spectroscopic pair, $`49.5\pm 3.3`$ km s-1. However, the observed difference in $`V_0`$ in HD 16429 Aa and Ab1 is probably due to differing wind velocities of the two stars. Radial velocities of photospheric lines are known to reflect the outward acceleration of the line formation region (Conti, Garmany, & Hutchings, 1977; Bohannan & Garmany, 1978). LS I +$`61^{}303`$ has a systemic velocity $`V_0=55\pm 6`$ km s-1 (Hutchings & Crampton, 1981). No radial velocities are available for HD 16429 C or ALS 7378. The proper motion of each star is from the Simbad database. For HD 16429 Aa and Ab1, $`V`$ and $`BV`$ were found using the values from Hiltner (1956) and the flux ratio determined above. For LS I +$`61^{}303`$, those values were found by averaging those by Drilling (1975) and Coyne & MacConnell (1983). LS I +$`61^{}303`$ has a somewhat redder $`BV`$ than the other stars, but it has previously been identified with the Cas OB6 association by Steele et al. (1998). $`V`$ and $`BV`$ for the remaining stars are from the Simbad database. Based on these similarities, HD 16429 A and its neighbors probably make up a small sub-cluster within the Cas OB6 association.

The stars in this sub-cluster were probably born in the same generation of star formation and therefore have approximately the same age and distance. Garmany & Stencel (1992) used main sequence fitting to determine a distance of $`2.40\pm 0.18`$ kpc to the Cas OB6 association. However, Frail & Hjellming (1991) used a kinematical model of the Galaxy to determine the distance to LS I +$`61^{}303`$ based on the velocities of H I absorption lines in its spectrum, and they placed that system at a closer distance of $`2.0\pm 0.2`$ kpc (confirmed by Steele et al. (1998)). Therefore I adopted a distance of 2.0 kpc for the group. Placing the stars on evolutionary tracks in a color-magnitude diagram is difficult since the reddening is variable across the association (Hanson & Clayton, 1993). For HD 16429 A, the spectral type of each component was used to obtain its intrinsic $`BV`$ (Wegner, 1994), which was then combined with the observed $`BV`$ for the system (Hiltner, 1956) to determine its reddening, $`E(BV)`$, and absolute magnitude, $`M_V`$. Howarth (1983) provides $`E(BV)=0.75`$ for LS I +$`61^{}303`$. The reddening of HD 16429 C and ALS 7378 was estimated assuming that each is a B0 V star (based on their $`B`$ and $`V`$ magnitudes relative to HD 16429 A) and following the same procedure used for HD 16429 A. Table 3 also lists the $`E(BV)`$ and $`M_V`$ for each star in the sub-cluster.

Böhm-Vitense (1981) provides temperatures of massive stars based on their spectral types, and I used her calibration to find their temperatures and bolometric corrections (Howarth & Prinja, 1989). Their luminosities were found using the expression log $`L/L_{\mathrm{}}=0.4(M_{bol}4.74)`$. The evolutionary tracks for non-rotating stars of Schaller et al. (1992), based on log $`T`$ and log $`L/L_{\mathrm{}}`$, provide estimates of the current and initial stellar masses as well as stellar ages, and these are also included in Table 3. It is interesting to note that the initial mass of HD 16429 Aa places a lower limit on the initial mass of the neutron star progenitor in LS I +$`61^{}303`$. It must have originated in a star with more than 39 $`M_{\mathrm{}}`$ when the group formed $`24`$ Myr ago to have evolved more quickly than HD 16429 Aa.

I especially thank Douglas Gies for his advice and comments regarding this manuscript. I also thank David Wingert and Wenjin Huang for taking many of the observations, and the KPNO staff for their assistance and travel support. I am also grateful to Dr. Alex Fullerton for sharing his period search code, and to Brian Mason and William Hartkopf of the U.S. Naval Observatory for sending me information from their speckle interferometric observations. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. Financial support was provided by the National Science Foundation through grant AST$``$0205297 (DRG). Institutional support has been provided from the GSU College of Arts and Sciences and from the Research Program Enhancement fund of the Board of Regents of the University System of Georgia, administered through the GSU Office of the Vice President for Research. I gratefully acknowledge all this support.