text:
# Cosmological Tests for a Linear Coasting Cosmology

## 1 Introduction

Standard cold dark matter (CDM) cosmology presents serious theoretical and observational problems as a model for an acceptable description of the Universe. This has motivated a search for alternative cosmological models . The first problem comes from the conflict between age of the Universe and the age of the oldest stars in Galactic globular clusters. The age constraints from old galaxies at high redshifts render “the age problem” even more acute .

Another major difficulty is the cosmological constant ($`\mathrm{\Lambda }`$) problem. While it is difficult to have a decent theoretical justification for $`\mathrm{\Lambda }`$ , there are serious doubts on the compatibility of the Hubble Deep Field (HDF) data with large $`\mathrm{\Lambda }`$. Initial analysis of Maoz and Rix  placed a somewhat stringent upper limit of $`\mathrm{\Omega }_\mathrm{\Lambda }<0.7`$ which have been somewhat softened by later analysis of Cooray et al. .

These problems have generated a lot of interest in an open FRW model with a linear evolution of the scale factor, $`a(t)t`$. In such a cosmology the universe expands with constant velocity; hence the term coasting cosmology . Notable among such models is a recent idea of Allen , in which such a scaling results in an $`SU(2)`$ cosmological instanton dominated universe. The Weyl gravity theory of Manheim and Kazanas  makes space for yet another possibility. Here again the FRW scale factor approaches a linear evolution at late times.

The need for investigating such a model comes from several considerations. Particle horizons occur in models with $`a(t)t^\alpha `$ for $`\alpha <1`$ and, therefore, linear coasting model does not suffer the horizon problem. Also, linear evolution of a scale factor is supported in alternative gravity theories (eg. non-minimally coupled scalar-tensor theories), where it turns out to be independent of the matter equation of state (see  and references therein). The scale factor in such theories does not constrain the matter density parameter and, therefore, does not present any flatness problem. Moreover, the age of the coasting universe is 1.5 times the age of standard CDM universe , which makes this model comfortably concordant with the ages of globular clusters. Finally, a linear coasting cosmology, independent of the equation of state of matter, is a generic feature in a class of models that attempt to dynamically solve the $`\mathrm{\Lambda }`$ problem . Such models have a scalar field non-minimally coupled to the curvature of the universe. With the evolution of time, the non-minimal coupling diverges, the scale factor quickly approaches linearity and the non-minimally coupled field acquires a stress energy that cancels the vacuum energy in the theory.

Interestingly it was noted by Perlmutter et al.  that the curve for $`\mathrm{\Omega }_\mathrm{M}=\mathrm{\Omega }_\mathrm{\Lambda }=0`$ (for which the scale factor would have linear evolution) is “practically identical to the best fit plot for an unconstrained cosmology”. Recently, it was shown by Dev et al.  that open linear coasting cosmology presents a good fit to the SNe Ia data. It was also demonstrated that this model is consistent with the primordial nucleosynthesis .

In this paper we consider constraints on the index $`\alpha `$ of a power law cosmology, $`a(t)t^\alpha `$, from two different tests: gravitational lensing statistics and age estimates of old high-redshift galaxies. The expected frequency of multiple imaging lensing events is a sensitive probe for the viability of a given cosmology. In view of the successful results of the above mentioned works , we used this test to constrain the power index $`\alpha `$ of the scale factor. Expected number of lens systems depends upon the index $`\alpha `$ through the angular diameter distances. By varying $`\alpha `$, the number of lenses changes which on comparison with the observations gives us the constrain on $`\alpha `$. Age measurements of old high-redshift galaxies give lower bound on the power index $`\alpha `$. This is based on the fact that the age of the universe in a given redshift is greater than or at least equal to the age of its objects. Since the age of the universe is a function of $`\alpha `$, we find the value of $`\alpha `$ which permits the existence of these old galaxies.

In Section 2 we introduce the ansatz for the power law cosmology and derive the angular diameter distance formula. In Section 3 we describe the tests and the constraints they present on the power index for that cosmology. We summarize our results and present the result of a joint constraint from various observational tests in Section 4.

## 2 Linear Coasting Cosmology

We consider a general power law cosmology with the scale factor given in terms of two arbitrary dimensionless parameters $`B`$ and $`\alpha `$

$$a(t)=B\frac{c}{H_0}\left(\frac{t}{t_0}\right)^\alpha ,$$
(1)
for an open FRW metric

$$ds^2=c^2dt^2a^2(t)\left[\frac{dr^2}{1+r^2}+r^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)\right].$$
(2)
Here $`t`$ is cosmic proper time and $`r,\theta ,\phi `$ are comoving spherical coordinates.

The expansion rate of the universe is described by a Hubble parameter, $`H(t)=\dot{a}/a=\alpha /t`$. The present expansion rate of the universe is defined by a Hubble constant, equal in our model to $`H_0=\alpha /t_0`$ (here and subsequently the subscript 0 on a parameter refers to its present value). The scale factor and the redshift are related to their present values by $`a/a_0=(t/t_0)^\alpha `$. As usual, the ratio of the scale factor at the emission and absorption of a null ray determines the cosmological redshift $`z`$ by

$$\frac{a_0}{a(z)}=1+z,$$
(3)
and the age of the universe is

$$t(z)=\frac{\alpha }{H_0(1+z)^{1/\alpha }}.$$
(4)
Using (3), we define the dimensionless Hubble parameter

$$h(z)\frac{H(z)}{H_0}=(1+z)^{1/\alpha }.$$
(5)
The present ‘radius’ of the universe is defined as (see Eq. 1)

$$a_0=B\frac{c}{H_0}.$$
(6)
In terms of the parameters $`\alpha `$ and $`B`$, the angular diameter distance between two different redshifts is:

$$D_\mathrm{A}(z_1,z_2,\alpha )=\frac{Bc}{(1+z_2)H_0}\mathrm{sinh}\left[\frac{1}{B}\frac{\alpha }{\alpha 1}\left\{(1+z_2)^{\frac{\alpha 1}{\alpha }}(1+z_1)^{\frac{\alpha 1}{\alpha }}\right\}\right].$$
(7)
In a limiting case, $`\alpha 1`$, we obtain the following expression

$$D_\mathrm{A}(z_1,z_2)=\frac{Bc}{2H_0}\frac{\left[(1+z_2)^{2/B}(1+z_1)^{2/B}\right]}{(1+z_1)^{1/B}(1+z_2)^{\frac{B+1}{B}}}.$$
(8)
The look-back time, which is the difference between the age of the universe when a particular light ray was emitted and the age of the universe now, we find as

$$\frac{cdt}{dz_\mathrm{L}}=\frac{c}{H_0(1+z_\mathrm{L})^{\frac{\alpha +1}{\alpha }}}.$$
(9)
## 3 Testing the model against observations

### 3.1 Gravitational lensing statistics

We consider a sample of 867 ($`z>1`$) high luminosity optical quasars which include 5 lensed quasars ($`1208+1011`$, H $`1413+117`$, LBQS $`1009+0252`$, PG $`1115+080`$, $`0142+100`$). This sample is taken from optical lens surveys such as the HST Snapshot survey , the Crampton survey , the Yee survey , Surdej survey , the NOT Survey  and the FKS survey . The lens surveys and quasar catalogs usually use V magnitudes, so we transform $`m_V`$ to a B-band magnitude using an average $`BV`$ colour of $`0.2`$ mag .

The differential probability $`d\tau `$ of a beam having a lensing event in traversing $`dz_\mathrm{L}`$ is

$$d\tau =n_0(1+z_L)^3\sigma \frac{cdt}{dz_\mathrm{L}}dz_\mathrm{L},$$
(10)
where $`n_0`$ is the present comoving number density of the lenses, $`\sigma `$ is the cross-section for lensing event and the quantity $`cdt/dz_\mathrm{L}`$ is given by (9).

For simplicity we use the Singular Isothermal Sphere (SIS) model for the lens mass distribution. The cross-section for lensing events for the SIS model is given by

$$\sigma =16\pi ^3\left(\frac{v}{c}\right)^4\left(\frac{D_{\mathrm{OL}}D_{\mathrm{LS}}}{D_{\mathrm{OS}}}\right)^2,$$
(11)
where $`v`$ is the velocity dispersion of the dark halo of the lensing galaxy. We define $`D_{\mathrm{OL}},D_{\mathrm{OS}}`$ and $`D_{\mathrm{LS}}`$ as the angular diameter distances from the observer to the lens, to the source and between the lens and the source, respectively.

Assuming no evolution of the galaxies, the comoving number density is modeled by Schechter function as

$$\mathrm{\Phi }(L,z=0)dL=\varphi _{}\left(\frac{L}{L_{}}\right)^{\stackrel{~}{\alpha }}\mathrm{exp}\left(\frac{L}{L_{}}\right)\frac{dL}{L_{}},$$
(12)
where $`\varphi _{}`$, $`\stackrel{~}{\alpha }`$ and $`L_{}`$ are the normalization factor, the index of faint-end slope, and the characteristic luminosity, respectively. We also assume that the velocity dispersion of dark matter halo $`v`$ is related to the luminosity $`L`$ by the Faber-Jackson relation for E/S0 galaxies $`v=v_{}\left(L/L_{}\right)^{1/\gamma }`$.

The differential probability is given by

$$d\tau =F^{}(1+z_\mathrm{L})^3\frac{H_0}{c}\left(\frac{H_0D_{\mathrm{OL}}D_{\mathrm{LS}}}{cD_{\mathrm{OS}}}\right)^2\frac{cdt}{dz_\mathrm{L}}dz_\mathrm{L},$$
(13)
where

$$F^{}=\frac{16\pi ^3}{cH_0^3}\varphi _{}v_{}^4\mathrm{\Gamma }\left(\stackrel{~}{\alpha }+\frac{4}{\gamma }+1\right)$$
(14)
is the dimensionless quantity, which measures the effectiveness of matter in producing multiple images. Table 1 lists Schechter and lens parameters for E/S0 galaxies, as suggested by  (hereafter $`LPEM`$ parameters). We neglect the contribution of spirals as lenses, as their velocity dispersion is small in comparison to E/S0 galaxies.

The constraints obtained on the cosmological paramters are highly dependent on the choice of lens and Schechter paramters. The lens and Schechter parameters should be determined in a highly correlated manner from a galaxy survey. The use of parameters derived from various surveys might introduce error. We consider LPEM parameters as they form one such set of paramters and they also take into account the morphological distribution of the E/S0 galaxies. Recently several galaxy surveys have come up with a much larger sample of galaxies. This has improved our knowledge of the galaxy luminosity function. However, these surveys donot classify the galaxies by their morphological type .

The differential optical depth of lensing in traversing $`dz_\mathrm{L}`$ with angular separation between $`\varphi `$ and $`\varphi +d\varphi `$ reads as

$`{\displaystyle \frac{d^2\tau }{dz_\mathrm{L}d\varphi }}d\varphi dz_\mathrm{L}=`$
$`F^{}(1+z_\mathrm{L})^3{\displaystyle \frac{H_0}{c}}{\displaystyle \frac{\gamma /2}{\mathrm{\Gamma }\left(\stackrel{~}{\alpha }+1+\frac{4}{\gamma }\right)}}{\displaystyle \frac{cdt}{dz_\mathrm{L}}}\left({\displaystyle \frac{H_0D_\mathrm{L}D_{\mathrm{LS}}}{cD_\mathrm{S}}}\right)^2\left({\displaystyle \frac{D_\mathrm{S}}{D_{\mathrm{LS}}}}\varphi \right)^{\frac{\gamma }{2}(\stackrel{~}{\alpha }+1+\frac{4}{\gamma })}\mathrm{exp}\left[\left({\displaystyle \frac{D_\mathrm{S}}{D_{\mathrm{LS}}}}\varphi \right)^{\frac{\gamma }{2}}\right]{\displaystyle \frac{d\varphi }{\varphi }}dz_\mathrm{L},`$
where $`\varphi =\mathrm{\Delta }\theta /8\pi (v_{}/c)^2`$ and $`v_{}`$ the velocity dispersion corresponding to the characteristic luminosity $`L_{}`$ in (12).

We make two corrections to the optical depth to get the lensing probability: magnification bias and selection function. Magnification bias, $`𝐁(m,z)`$, is to take into account the increase in the apparent brightness of a quasar due to lensing, which, in turn, increases the expected number of lenses in flux limited sample.

The bias factor for a quasar at redshift $`z`$ with apparent magnitude $`m`$ is given by

$$𝐁(m,z)=M_0^2\text{B}(m,z,M_0,M_2),$$
(16)
where

$$\text{B}(m,z,M_1,M_2)=2\left(\frac{dN_\mathrm{Q}}{dm}\right)^1_{M_1}^{M_2}\frac{dM}{M^3}\frac{dN_\mathrm{Q}}{dm}(m+2.5\mathrm{log}(M),z).$$
(17)
In the above equation $`\left(dN_\mathrm{Q}(m,z)/dm\right)`$ is the measure of number of quasars with magnitudes in the interval $`(m,m+dm)`$ at redshift $`z`$. We can allow the upper magnification cutoff $`M_2`$ to be infinite, though in practice we set it to be $`M_2=10^4`$. $`M_0`$ is the minimum magnification of a multiply imaged source and for the SIS model $`M_0=2`$.

We use Kochanek’s “best model” (K96) for the quasar luminosity function:

$$\frac{dN_\mathrm{Q}}{dm}(m,z)\left(10^{a(m\overline{m})}+10^{b(m\overline{m})}\right)^1,$$
(18)
where the bright-end slope $`a`$ and faint-end slope $`b`$ are constants, and the break magnitude $`\overline{m}`$ evolves with redshift:

$$\overline{m}=\{\begin{array}{cc}m_o+(z1)\hfill & \text{for }z<1,\hfill \\ m_o\hfill & \text{for }1<z3,\hfill \\ m_o0.7(z3)\hfill & \text{for }z>3.\hfill \end{array}$$
(19)
Fitting this model to the quasar luminosity function data in  for $`z>1`$, Kochanek finds that “the best model” has $`a=1.07\pm 0.07`$, $`b=0.27\pm 0.07`$ and $`m_o=18.92\pm 0.16`$ at B magnitude. The magnitude corrected probability, $`p_i`$, for the quasar $`i`$ with apparent magnitude $`m_i`$ and redshift $`z_i`$ to get lensed is:

$$p_i=\tau (z_i)𝐁(m_i,z_i).$$
(20)
Selection effects are caused by limitations on dynamic range, limitations on resolution and presence of confusing sources such as stars. The survey can only detect lenses with magnifications larger than $`M_f`$. This sets the lower limit on the magnification. Therefore the $`M_1`$ in the bias function (17) gets replaced by $`M_f(\theta )`$ (for details, see K93), which is given as

$$M_f=M_0(f+1)/(f1),$$
(21)
with

$$f=10^{0.4\mathrm{\Delta }m(\theta )}.$$
(22)
The corrected lensing probability and image separation distribution function for a single source at redshift $`z_\mathrm{S}`$ are given in K96

$$p_i^{^{}}(m,z)=p_i\frac{d(\mathrm{\Delta }\theta )p_c(\mathrm{\Delta }\theta )\text{B}(m,z,M_f(\mathrm{\Delta }\theta ),M_2)}{\text{B}(m,z,M_0,M_2)},$$
(23)
and

$$p_{ci}^{^{}}=p_{ci}(\mathrm{\Delta }\theta )\frac{p_i}{p_i^{^{}}}\frac{\text{B}(m,z,M_f(\mathrm{\Delta }\theta ),M_2)}{\text{B}(m,z,M_0,M_2)},$$
(24)
where

$$p_c(\mathrm{\Delta }\theta )=\frac{1}{\tau (z_S)}_0^{z_S}\frac{d^2\tau }{dz_Ld(\mathrm{\Delta }\theta )}𝑑z_L.$$
(25)
Equation (24) defines the configuration probability. It is the probability that the lensed quasar $`i`$ is lensed with the observed image separation.

To get selection function corrected probabilities, we divide our sample into two parts—the ground based surveys and the HST survey. We use the selection functions as suggested in K93.

In our present calculations we do not consider the extinction effects due to the presence of dust in the lensing galaxies.

Finally, we use the above basic equations to perform the following tests :

* (i) The sum of the lensing probabilities $`p_i^{^{}}`$ for the optical QSOs gives the expected number of lensed quasars, $`n_\mathrm{L}=p_i^{^{}}`$. The summation is over the given quasar sample.We look for those values of the parameter for which the adopted optical sample has exactly five lensed quasars (that is those values of the parameters for which $`n_\mathrm{L}=5`$).

We start with a two parameter fit. We allow $`\alpha `$ to vary in the range ($`0.0\alpha 2.0`$) and $`B`$ to vary in the range ($`0.5B10.0`$). We observe that for $`B1`$, $`D_A`$ becomes independent of it. It can be easily checked from equation (7) that for large values of $`B`$, $`\mathrm{sinh}(x)x`$ making $`D_A`$ independent of $`B`$. From the previous works constraining power law cosmology , the value of $`B=1`$ is found to be compatible with observations. Incidentally, we can estimate the present scale factor of the universe as $`a_0c/H_0`$, hence we use $`B=1`$ in further analysis.

Fig. 1 shows the predicted number of lensed quasars for the above specified range of $`\alpha `$. We obtain $`n_\mathrm{L}=5`$ for $`\alpha =1.06`$. We further generate $`10^4`$ quasar samples (each sample has 867 quasars) using bootstrap method and find best fit $`\alpha `$ for each data set to obtain error bars on $`\alpha `$. We finally obtain $`\alpha =1.09\pm 0.3`$.
* (ii) We also perform maximum likelihood analysis to determine the value of $`\alpha `$, for which the observed sample becomes the most probable observation. The likelihood function is

$$=\underset{i=1}{\overset{N_\mathrm{U}}{}}(1p_i^{^{}})\underset{k=1}{\overset{N_\mathrm{L}}{}}p_k^{^{}}p_{ck}^{^{}}.$$
(26)
Here $`N_\mathrm{L}`$ is the number of multiple-imaged lensed quasars, $`N_\mathrm{U}`$ is the number of unlensed quasars, $`p_k^{^{}}`$, the probability of quasar $`k`$ to get lensed is given by Eq. 23 and $`p_{ck}^i`$, the configuration probability, is given by Eq. 24. The best fit ($`_{\mathrm{max}}`$) occurs for $`\alpha =1.13`$. We see that $`0.85\alpha 1.56`$ at $`1\sigma `$ ($`68\%`$ confidence level) and $`0.65\alpha 2.33`$ at $`2\sigma `$ ($`95.4\%`$ confidence level).

However, gravitational lensing statistics is susceptible to a number of uncertainties. The constraints obtained on the cosmological parameters from the statistics of strong lensing may vary after inclusion of the uncertainties in the luminosity function, lensing cross-section for galaxies (E/S0) and quasar luminosity function, role of spirals and the dust extinction.

### 3.2 Constraints from age estimates of high-z galaxies

Another observational test which can constraint $`\alpha `$ is the age measurement of the old high redshift galaxies (OHRG). These constraints are more stringent than those obtained from globular cluster age measurements . Here, we consider the galaxy 3C65 at $`z=1.175`$ (4 Gyr old) , at $`z=1.55`$ (3.5 Gyr old; 53W091)  and a 4 Gyr old galaxy 53W069 at $`z=1.43`$ . These are the minimum ages of these galaxies as indicated by best fitting spectral synthesis models. The age of the universe at a given redshift is greater than or at least equal to the age of its oldest objects at that redshift. In power law cosmologies the age of the universe increases with increasing $`\alpha `$. Hence this test gives lower bound on $`\alpha `$. This can be checked if we define the dimensionless ratio:

$$\frac{t(z)}{t_\mathrm{g}}=\frac{f(\alpha ,z)}{H_0t_\mathrm{g}}1.$$
(27)
The $`t_\mathrm{g}`$ is the age of an old object and $`f(\alpha ,z)=\alpha /(1+z)^{1/\alpha }`$. The error bar on $`H_0`$ determines the extreme value of $`t_\mathrm{g}`$. The lower limit on $`H_0`$ was recently updated to nearly $`10\%`$ of accuracy by Freedman : $`H_0=70\pm 7`$ km/sec/Mpc; $`1\sigma `$. The constraints from SNe Ia data also point to $`H_0>60`$ km/sec/Mpc . For the galaxy 3C65, the lower limit on age (4.0 Gyr) yields: $`0.26H_0t_\mathrm{g}0.32`$. Similarly for the galaxy 53W069 we have $`0.26H_0t_\mathrm{g}0.32`$. For the galaxy 53W091 at $`z=1.55`$ the lower limit on age (3.5 Gyr) gives $`0.23H_0t_\mathrm{g}0.28`$.

Fig. 2 shows the variation of the function $`f(\alpha ,z)`$ with the redshift $`z`$ for several values of $`\alpha `$. We see that $`\alpha `$ should be at least 0.8, in order to allow for these OHRG to exist.

## 4 Summary and a combined constraint

The main results of the present paper along the with constraints obtained from the SNe Ia data  are summarized in Table 2. The motivation for our work was to establish the viability of a linear coasting cosmology $`a(t)=t`$. Using gravitational lensing statistics, we find that such a coasting is accommodated within $`1\sigma `$. The age determination of OHRG gives as a lower bound $`\alpha 0.8`$ for a power law cosmology $`a(t)t^\alpha `$. Moreover, for such a cosmology $`H_0t_0=\alpha `$. With updated value of $`H_0=70\pm 7`$ km/sec/Mpc and $`t_0=14\pm 2`$ Gyr , it gives $`\alpha =0.98\pm 0.25`$. Dev et al (2001) reported that $`\alpha =1.0`$ is consistent with SNe Ia data (within 68% confidence level). We find that $`\alpha =1.0`$ is in concordance with the listed observational tests. It is interesting to observe that the lensing analysis barely accomodates an Einstein-de Sitter universe ($`\alpha =2/3`$) at $`2\sigma `$. Similarly, the age determination of OHRG also rules out $`\alpha =2/3`$. We conclude that the coasting cosmology with strictly linear evolution of scale factor, $`a(t)=t`$, cannot be ruled out on the basis of these observations.

## Acknowledgments

The authors are grateful to T. D. Saini, A. Habib and N. Mahajan for useful discussions during the course of this work. We also thank the referee for the useful comments.


format: text

text:
# Spatial resolution of spin waves in an ultra-cold gas

## Abstract

We present the first spatially resolved images of spin waves in a gas. The complete longitudinal and transverse spin field as a function of time and space is reconstructed. Frequencies and damping rates for a standing-wave mode are extracted and compared with theory.

Macroscopic collective behavior can arise in a quantum gas even above the onset of degeneracy. When indistinguishable atoms collide, the scattering can be significantly altered by the need to symmetrize the scattered waves. Examples of quantum scattering effects are seen in the strong polarization dependence of heat conduction and viscosity coefficients in spin polarized 3He emery1964  and of thermalization rates in nondegenerate Fermi gases jin1999 . Quantum collisional effects can also produce spatio-temporal spin oscillations known as spin waves. In a recent experiment, we observed in an ultra-cold atomic gas a startling spin migration effect which we called “anomalous spin-state segregation” lewando2002 . Several theory groups have shown that this effect can be thought of as a half-cycle of a large-amplitude, over-damped spin wave newswtheory . In this Letter we verify experimentally their viewpoint and present spatially resolved images of coupled transverse and longitudinal spin perturbances propagating through a magnetically trapped atomic cloud. Frequencies and damping rates for the standing-wave excitations are studied as a function of density and temperature.

Although the concept of spin waves in ferromagnets dates back to Bloch’s predictions in 1930 bloch1930 , the theory of spin waves in liquids and in dilute gases was not formulated until the much later spinwavetheory . The first evidence for spin waves (in the form of extra resonances in nuclear magnetic resonance spectra) occurred in spin polarized hydrogen and soon after in polarized 3He and dilute mixtures of 3He in 4He spinwaveexp . On a microscopic level, spin waves arise from interference effects in lowest partial-wave collisions. When identical particles collide, the $`180^{}`$ backscattering event causes a scattered atom to propagate along a trajectory indistinguishable from that of either a forward scattered atom or an unscattered atom. When the two atoms are spin aligned, the backward scattering contribution is indistinguishable from the forward scattering contribution, which adds a factor of two to the mean-field collisional energy shift. When the spins are antiparallel, the backward scattering event is distinguishable and the mean-field shift arises from the forward scattering event only. In the case of intermediate spin alignment, the backward scattering event is only partially distinguishable from the unscattered event, and one needs to add coherently two unparallel spin amplitudes, leading to a rotation of the spins of the scattered atoms. The cumulative effect of many such spin-rotating collisions is such that inhomogeneities of spin propagate like waves rather than diffusively. The characteristic frequency scale for the spin rotation effect is the exchange collision rate, $`\omega _{\text{exch}}=4\pi \mathrm{}an/m`$, where $`m`$, $`n`$, and $`a`$ are the mass of the atoms, the number density, and the s-wave scattering length respectively.

The experimental apparatus and spectroscopy method used to prepare and probe the ultra-cold gas are described in lewando2002  and briefly will be summarized here. 87Rb atoms are precooled and trapped in a magneto-optical trap, transferred to a Ioffe-Pritchard magnetic trap via a servo-controlled linear track, and cooled further using forced radio-frequency evaporation. The trapping frequencies are {7, 230, 230} Hz, and typical atom cloud parameters are $`n10^{13}`$ cm-3 and $`T600`$ nK, several times the Bose-Einstein condensation temperature. The radial frequency is high enough that for the effects considered here it is sufficient to average the cloud over the radial dimensions, effectively reducing the cloud to a one dimensional density distribution. There are two conveniently trapped atomic hyperfine sates, labelled $`|1`$ and $`|2`$, which together make up a pseudo-spin doublet states . The spin state is manipulated using a two-photon coupling transition lewando2002 . Typical Rabi frequencies are around 3 kHz.

The transition frequency $`\mathrm{\Delta }(z)`$ is in general a function of the axial position $`z`$, due to the spatial inhomogeneity of the Zeeman shift and of the density-dependent mean-field shift. In Ref. lewando2002  we describe in detail how the bias magnetic field can be tuned so as to cause the Zeeman shifts to cancel in whole or in part the mean-field shifts. As a consequence, the spatial curvature of the frequency, $`^2\mathrm{\Delta }(z)/z^2`$, is a controllable parameter of the experiment which may be tuned in real time, even as a wave is propagating.

To describe the effective two-level system, we use the language of the Bloch sphere (see allen1987 ) and take the axis *w* as the “longitudinal” population inversion, and the axes *u* and *v* as the “transverse” coherences sxsz . In a typical experiment, the atom-preparation cycle concludes with the atoms initially all in the $`|1`$ state, along the *w*-axis. Applying a $`\pi /2`$ pulse starts the evolution of the wave by rotating all the spins to lie along the *v*-axis. When $`\mathrm{\Delta }(z)`$ is set to have nonzero curvature, the spins along the long axis of the cloud begin to fan out in the *u-v* plane as the relative phase of the $`(|1+|2)`$ coherent superposition develops a spatial dependence. Optionally at this point in the evolution, we can tune $`^2\mathrm{\Delta }(z)/z^2`$ back to zero or leave it unchanged and watch the spins evolve in an inhomogeneous potential. In either case, the atoms’ thermal motion now carries spin information back and forth along the length of the cloud, and spin-rotating collisions can tilt the spins up out of the *u-v* plane meanfield . Eventually, spin-currents develop that drive oscillations about an equilibrium spin distribution. Our driving potential, determined by the curvature of $`\mathrm{\Delta }(z)`$, is even in $`z`$, and we predominantly drive the lowest order symmetric mode: a two-node standing wave.

The longitudinal and transverse components of the Bloch vectors are measured as a function of time and space. The method for projecting the longitudinal component is similar to that described in lewando2002 . After an initial resonant $`\pi /2`$ pulse, the Bloch vectors are allowed to evolve for times up to $`t=800`$ ms, and then the atoms are imaged, thereby projecting the Bloch vectors along the cloud onto the two spin states $`|1`$ and $`|2`$. The cloud is radially averaged and broken up into 25 equal bins along the axial direction for spatial resolution of the spin dynamics. The $`|1`$ and $`|2`$ populations are measured on separate experimental shots, and the longitudinal angle $`\theta `$ of the Bloch vector is extracted. ($`\theta =\pi /2`$ corresponds to equal populations, and $`\theta =0`$ or $`\pi `$ correspond to all the atoms in $`|1`$ or $`|2`$ respectively.) The evolution is traced out by repeating the experiment for many values of evolution time. Fig. 1a shows the time evolution of the longitudinal component of the spin across the cloud. The initial oscillation appears to be “spin-state segregation” lewando2002  but in fact is the first half-cycle of a spin wave as the Bloch vectors in the cloud center tip up from the transverse plane and the vectors near the cloud edge tip down.

The transverse spin measurement consists of a two-pulse $`\pi /2\pi /2`$ Ramsey pulse sequence separated by an evolution time $`t`$ ramsey1956 . For each time evolution time $`t`$ and for each spatial bin along the trap axis, we measure the final density of the atoms in $`|1`$, resulting in a set of roughly sinusoidal (but with slightly time-varying phase and amplitude) Ramsey fringes for each location in space. The transverse phase and amplitude of the Bloch vector are extracted from the fringes at each time $`t`$ for each radially-averaged section of the cloud by performing a sinusoidal fit over a small window centered about $`t`$ in each spatial bin transverse . Fig. 1b shows typical phases extracted with this method. The spatial and temporal oscillations around the mean phase are the transverse projections of the spin waves.

The complete trajectories of the Bloch vector can be determined for each spatial bin by combining the transverse and longitudinal components. Normalizing the length of each vector to the length at time $`t=0`$ removes any effects of decoherence (*e.g.* from inelastic loss and magnetic field fluctuations) and allows for a study of the equilibrium phases of the Bloch vectors as a function of the position across the cloud for various $`\mathrm{\Delta }(z)`$. Fig. 2 shows the trajectories of the reconstructed Bloch vectors in each spatial bin. Each bin has been offset by 1 rad in the transverse direction from the adjacent bins for clarity and to depict schematically the position across the cloud. The oscillations can be fit to extract spin wave frequency, amplitude, and damping rate.

The steady-state Bloch vector arrangement is determined by the curvature of the potential. In order to equilibrate, the Bloch vectors across the cloud must arrange themselves so that the torque from the driving potential $`\mathrm{\Delta }(z)`$ (which tends to spread the vectors in the *u-v* plane) is balanced by a torque from the local curvature in the spin field, which is spread out in the *v-w* plane. In the presence of a nonuniform potential, the equilibrium configuration of the Bloch vectors is curved longitudinally across the cloud to match the curvature of $`\mathrm{\Delta }(z)`$, *i.e.* the Bloch vectors lie in the *v-w* plane but are not all collinear with the *v*-axis. However, if the bias magnetic field is chosen to cancel the density shift, then $`\mathrm{\Delta }(z)`$ is predominantly flat across the cloud, and the equilibrium corresponds to the Bloch vectors everywhere in the cloud collinear with the *v*-axis. These equilibrium configurations are illustrated in Fig. 2a and Fig. 2b, which show spin-wave evolution in curved and flat potentials respectively.

Next, the density dependence of the excitation frequency was studied using transverse spin waves. The excitation frequency exhibits a strong density dependence, scaling roughly as $`1/n`$, shown in Fig. 3a. The linear regime for these excitations occurs only at extremely small amplitude, where spin waves are difficult to observe due to signal-to-noise limitations. In the nonlinear regime there is a dependence of frequency on amplitude. In order to remove this dependence, we excite spin waves over a range of amplitudes by controlling the curvature of the potential, $`^2\mathrm{\Delta }(z)/z^2`$, and we extrapolate the frequency to the low amplitude limit for each density. The solid line is a numerical calculation predicting the frequency obtained from solving the one-dimensional Boltzmann transport equations williams2002b . The two dotted lines represent the weak and strong coupling limits, $`\omega _{\text{exch}}\omega _{}`$ and $`\omega _{\text{exch}}\omega _{}`$ respectively:

$`\omega _\text{w}`$ $`=`$ $`2\omega _{}{\displaystyle \frac{\omega _{\text{exch}}}{2}}`$ (1)
$`\omega _\text{s}`$ $`=`$ $`{\displaystyle \frac{k_{\text{eff}}^2(k_\text{b}T)}{m\omega _{\text{exch}}}}.`$ (2)
where $`\omega _{}`$ is the axial trap frequency and $`k_{\text{eff}}`$ is an effective wave vector. The spatial wave vector $`k=2\pi /\lambda `$ of the spin waves is studied by finding the distance $`\lambda /2`$ between the spatial nodes, *i.e.* the bins where the orientation of the Bloch vector remains constant. The wave vectors $`k`$ and $`k_{\text{eff}}`$ differ by a numerical factor of order unity due to the finite size of the sample. For all the work discussed here, the lowest order symmetric mode, $`\eta =2`$, is excited by the symmetric inhomogeneous potential mode . The wavelength of the $`\eta =2`$ mode is observed to be directly proportional to the size of the trap. As predicted by Eqs. 1 and 2, there is no temperature dependence for the excitation frequency, since we observe $`k`$ to scale inversely as the size of the cloud, *i.e.* as $`T^{1/2}`$. The numerical calculation agrees well with experiment (see Fig. 3a).

The Bloch vectors initially precess about the equilibrium value with some initial amplitude. As the Bloch vectors precess, their trajectories spiral towards the equilibrium configuration due to damping of the spin wave excitation. The damping is predicted to be a combination of two effects: spin diffusion and Landau damping levitov1999 . Spin diffusion dominates in the high density, or hydrodynamic, regime. For intermediate densities, there is a maximum in the damping rate due to the combination of the two effects. It near this intermediate regime that these experiments are performed. We observe a relationship between spin wave lifetime and density which is good agreement with theory (see Fig. 3b).

Our original motivation to explore spin dynamics in uncondensed clouds arose from our plan to extend to finite temperatures our earlier studies of spin coherence in a condensate hall1998b  and in particular to understand the extent to which an unpolarized normal cloud acts as a decohering thermal reservoir. The existence of spin waves will complicate these studies. A pure condensate cannot support spin waves, but it will trivially separate due to small differences in the scattering lengths hall1998 . Spin oscillations in the thermal cloud can have a considerable effect on the spin dynamics of a finite temperature condensate, and conversely, having a small stationary cluster of condensed atoms will significantly perturb spin waves in a thermal cloud. Such a two component system is a rich system for studying spin dynamics in ultra-cold gases levitov1999 .

In conclusion, we have observed collective spin excitations in a nondegenerate ultra-cold gas. Particle indistinguishability and the need to symmetrize collisions between like particles in coherent superpositions give rise to spatial and temporal oscillations of the Bloch vectors across a cloud of trapped atoms. We have separately imaged both transverse and longitudinal components of these spin waves and studied the effects of temperature and density on the spin wave frequency, paving the way for studying spin dynamics in a finite-temperature condensate system.

We acknowledge useful conversations with the other members of the JILA BEC collaboration. This work was supported by grants from the NSF and NIST.


text:
# Partonic Effects on Pion Interferometry at the Relativistic Heavy Ion Collider

## Abstract

Using a multiphase transport (AMPT) model that includes both initial partonic and final hadronic interactions, we study the pion interferometry at the Relativistic Heavy Ion Collider. We find that the two-pion correlation function is sensitive to the magnitude of the parton scattering cross section, which controls the parton density at which the transition from the partonic to hadronic matter occurs. Also, the emission source of pions is non-Gaussian, leading to source radii that can be more than twice larger than the radius parameters extracted from a Gaussian fit to the correlation function.

Particle interferometry based on the Hanbury-Brown Twiss (HBT) effect has long been used to measure the size of an emission source hbt . In heavy ion collisions, it has been suggested that the HBT study can provide information not only on the spatial extent of the emission source but also on its expansion velocity and emission duration Pratt:su ; Bertsch:1988db ; Pratt:zq ; Rischke:1996em . In particular, the long emission time as a result of the phase transition from the quark-gluon plasma to hadronic matter in relativistic heavy ion collisions is expected to lead to an emission source which has a much larger radius in the direction of the total transverse momentum of detected two particles ($`R_{\mathrm{out}}`$) than that perpendicular to both this direction and the beam direction ($`R_{\mathrm{side}}`$) Rischke:1996em ; Soff:2000eh ; Soff:2001hc . Since the quark-gluon plasma is expected to be formed in heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC), it is thus surprising to find that the extracted ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ from a Gaussian fit to the measured two-pion correlation function in Au+Au collisions is close to one Adler:2001zd ; Johnson:2001zi ; Adcox:2002uc . Also, the extracted radius parameters are small compared to theoretical predictions based on the hydrodynamical model Soff:2000eh . These experimental results, especially the small value of $`R_{\mathrm{out}}/R_{\mathrm{side}}`$, have been attributed to strong space-time and momentum correlations in the emission source Tomasik:1998qt .

In this Letter, a multiphase transport model (AMPT) Zhang:2000bd ; Lin:2001cx ; Lin:2001zk , that includes both initial partonic and final hadronic interactions, is used to study the pion interferometry in central Au+Au collisions at RHIC at $`\sqrt{s}=130A`$ GeV. The AMPT model is a hybrid model that uses the minijet partons from hard processes and the strings from soft processes in the HIJING model Wang:1991ht  as the initial conditions for modeling the collision dynamics. The time evolution of partons is then modeled by the ZPC Zhang:1997ej  parton cascade model. At present, this model includes only parton-parton elastic scatterings with an in-medium cross section given by:

$`{\displaystyle \frac{d\sigma _p}{d\widehat{t}}}={\displaystyle \frac{9\pi \alpha _s^2}{2}}\left(1+{\displaystyle \frac{\mu ^2}{\widehat{s}}}\right){\displaystyle \frac{1}{(\widehat{t}\mu ^2)^2}},`$ (1)
where the strong coupling constant $`\alpha _s`$ is taken to be 0.47, and $`\widehat{s}`$ and $`\widehat{t}`$ are the Mandelstam variables. The effective screening mass $`\mu `$ depends on the temperature and density of the partonic matter but is taken as a parameter for fixing the magnitude and angular distribution of the parton scattering cross section. After these minijet partons stop interacting, they are combined with their parent strings to fragment to hadrons using the Lund string fragmentation model as implemented in the PYTHIA routine Sjostrand:1994yb . The final-state hadronic scatterings are then modeled by the ART model Li:1995pr .

The default AMPT model has been quite successful in describing the measured rapidity distributions of charge particles, particle to antiparticle ratios, and the spectra of low transverse momentum pions and kaons Lin:2001cx  in heavy ion collisions at the Super Proton Synchrotron (SPS) and RHIC. Since the initial energy density in Au+Au collisions at RHIC is expected to be much larger than the critical energy density at which the hadronic matter to quark-gluon plasma transition would occur Kharzeev:2001ph ; Zhang:2000nc , the AMPT model has been extended to convert the initial excited strings into partons Lin:2001zk . In this string melting scenario, hadrons, that would have been produced from string fragmentation, are converted instead to valence quarks and/or antiquarks. Interactions among these partons are again described by the ZPC parton cascade model. The transition of the partonic matter to the hadronic matter is, however, achieved using a simple coalescence model, which combines two nearest partons into mesons and three nearest partons into baryons or anti-baryons. Using parton scattering cross sections of 6-10 mb, the extended AMPT model is able to reproduce both the centrality and transverse momentum (below 2 GeV$`/c`$) dependence of the elliptic flow measured in Au+Au collisions at $`\sqrt{s}=130A`$ GeV at RHIC Ackermann:2000tr .

From the AMPT model, the source of emitted particles is obtained from their space-time coordinate $`x`$ and momentum $`𝐩`$ at freezeout, i.e., at their last interactions. Denoting the single-particle emission function for pions by $`S(x,𝐩)`$, the HBT correlation function for two identical pions in the absence of final-state interactions, such as the Coulomb interaction, is then given by Pratt:su ; Wiedemann:1999qn :

$`C_2(𝐐,𝐊)=1+{\displaystyle \frac{d^4x_1d^4x_2S(x_1,𝐊)S(x_2,𝐊)\mathrm{cos}\left[Q(x_1x_2)\right]}{d^4x_1S(x_1,𝐩_\mathrm{𝟏})d^4x_2S(x_2,𝐩_\mathrm{𝟐})}},`$ (2)
where $`𝐩_\mathrm{𝟏}`$ and $`𝐩_\mathrm{𝟐}`$ represent the momenta of the two pions, respectively, $`𝐊=(𝐩_\mathrm{𝟏}+𝐩_\mathrm{𝟐})/2`$, and $`Q=(𝐩_\mathrm{𝟏}𝐩_\mathrm{𝟐},E_1E_2)`$. Expecting that the emission function is sufficiently smooth in the momentum space, one can evaluate the correlation function by using $`𝐩_\mathrm{𝟏}`$ and $`𝐩_\mathrm{𝟐}`$ for $`𝐊`$ in the numerator of Eq.(2). Furthermore, the three-dimensional correlation function in $`𝐐`$ is usually shown as a function of the invariant relative momentum ($`Q_{\mathrm{inv}}=\sqrt{Q^2}`$) or as a function of the projection of $`𝐐`$ in the “out-side-long” ($`\mathrm{𝑜𝑠𝑙}`$) system Bertsch:1988db ; Pratt:zq , defined by the beam direction ($`Q_{\mathrm{long}}`$), the direction along the total momentum of the two particles in the transverse plane ($`Q_{\mathrm{out}}`$), and the direction orthogonal to the above two directions ($`Q_{\mathrm{side}}`$).

Using the emission function obtained from the AMPT model for central ($`b=0`$ fm) Au+Au collisions at $`\sqrt{s}=130A`$ GeV, we have evaluated the correlation function $`C_2(𝐐,𝐊)`$ in the longitudinally comoving system frame using the program Correlation After Burner pratt:uf . In Fig. 1, we show the one-dimensional projections of the calculated correlation function including final-state Coulomb interactions for midrapidity ($`0.5<y<0.5`$) charged pions with transverse momentum $`125<p_\mathrm{T}<225`$ MeV$`/c`$. Also shown are the measured $`\pi ^{}`$ correlation functions from central collisions by the STAR collaboration without removing the effect due to the Coulomb interaction Adler:2001zd . In evaluating the one-dimensional projections of the correlation function onto one of the $`Q_{\mathrm{out}},Q_{\mathrm{side}}`$ and $`Q_{\mathrm{long}}`$ axes, we have integrated the other two $`𝐐`$ components over the range $`035`$ MeV/$`c`$. The dash-dotted curves in Fig. 1 are results from the default AMPT model (no string melting) with a parton scattering cross section of $`\sigma _p=3`$ mb, while other curves are those from the extended AMPT model with string melting but different values for $`\sigma _p`$. It is seen that with string melting both the width of the $`Q_{\mathrm{inv}}`$ correlation function and its height decrease with increasing $`\sigma _p`$. The decreasing width can be more clearly seen from the calculated correlation functions without including final-state Coulomb interactions long , as the value of the correlation function in this case is exactly two at $`Q_{\mathrm{inv}}=0`$. Since the integration of the other two $`𝐐`$ components over a fixed range ($`035`$ MeV/$`c`$) gives a smaller value when the correlation function becomes narrower in $`𝐐`$, the height of the one-dimensional projection of the correlation function decreases with increasing $`\sigma _p`$. To reproduce the measured one-dimensional correlation functions by the STAR collaboration, we find that a parton scattering cross section of about 10 mb is required in the extended AMPT model with string melting. Since in our model partons are converted to hadrons after they make their last scattering, the parton scattering cross section thus controls the density at which partons are converted to hadrons. A larger $`\sigma _p`$ leads to a lower density and larger volume for the parton to hadron phase transition. We expect that the two-pion correlation data could also be reproduced by using a smaller parton scattering cross section but delaying the hadronization until the critical density is similar to the one given by the parton freezeout condition using a scattering cross section of 10 mb in the parton cascade model long .

The size of the emission source can be determined from the emission function via the curvature of the correlation function at $`𝐐=0`$:

$`R_{ij}(K)^2={\displaystyle \frac{1}{2}}{\displaystyle \frac{^2C_2(𝐐,𝐊)}{Q_iQ_j}}|_{𝐐=0}`$
$`=D_{x_i,x_j}(K)D_{x_i,\beta _jt}(K)D_{\beta _it,x_j}(K)+D_{\beta _it,\beta _jt}(K),`$ (3)
where $`x_i(i=13)`$ denotes the projections of the particle position at freezeout in the osl system, i.e., $`x_{\mathrm{out}}`$, $`x_{\mathrm{side}}`$ and $`x_{\mathrm{long}}`$, and $`\beta =𝐊/K_0`$ with $`K_0`$ being the average energy of the two particles. In the second line of Eq.(3), we have the variance $`D_{x,y}=xyxy`$, with $`x`$ denoting the average value of $`x`$.

Since only the correlation function is measured in experiments, the size of emission source is usually estimated by fitting the measured correlation function $`C_2(𝐐,𝐊)`$, that has been corrected for effects due to final-state Coulomb interactions, with a four-parameter Gaussian function:

$`C_2(𝐐,𝐊)=1+\lambda \mathrm{exp}\left({\displaystyle \underset{i=1}{\overset{3}{}}}R_{ii}^2(K)Q_i^2\right).`$ (4)
If the emission source is Gaussian in space-time, then for central heavy ion collisions considered here the radius parameters obtained from the above Gaussian fit to the correlation function would be the same as those determined directly from the emission function of the source via Eq.(3). However, because of space-time correlations in the emission function, such as those induced by the collective flow, the radius parameters in general do not give the true source size.

In Fig. 2, we show by solid curves the transverse mass ($`m_\mathrm{T}`$) dependence of the source radii $`R_{\mathrm{out}}`$ (upper-left panel), $`R_{\mathrm{side}}`$ (upper-right panel), and $`R_{\mathrm{long}}`$ (lower-left panel) determined from the emission function for midrapidity charged pions given by the AMPT model with string melting and parton cross section of $`\sigma _p=10`$ mb for central Au+Au collisions at $`\sqrt{s}=130A`$ GeV. The radius parameters determined from the Gaussian fit to the three-dimensional correlation function in $`𝐐`$ without including final-state Coulomb interactions are shown by curves with squares. For the $`m_\mathrm{T}`$ values shown here, these radius parameters are about a factor of 2 to 3 smaller than the source radii obtained directly from the emission function. The emission source from the AMPT model thus deviates appreciably from a Gaussian one. The radius parameters from the Gaussian fit to the correlation function have similar values as the experimental ones from a Gaussian fit to the measured correlation function after correcting for the final-state Coulomb interactions Adler:2001zd . The parameter $`\lambda `$ (scaled up by a factor of 10) from the Gaussian fit to the correlation function is shown in the lower-right panel by the curve with squares. It has a value of about 0.5 at low $`m_\mathrm{T}`$ but increases to about 1 at large $`m_\mathrm{T}`$. Both the source radii and fitted radius parameters decrease with increasing $`m_\mathrm{T}`$. On the other hand, these radii increase with increasing parton scattering cross section as a result of both a larger source size and stronger collective expansion long .

Since Eq.(3) gives

$`R_{\mathrm{out}}^2`$ $`=`$ $`D_{x_{\mathrm{out}},x_{\mathrm{out}}}2D_{x_{\mathrm{out}},\beta _{}t}+D_{\beta _{}t,\beta _{}t},`$ (5)
and $`R_{\mathrm{side}}^2=D_{x_{\mathrm{side}},x_{\mathrm{side}}}`$, the ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ contains information about the duration of emission and has been studied extensively Soff:2000eh ; Adler:2001zd . This connection between $`R_{\mathrm{out}}`$ and the emission duration becomes clearer if we neglect the $`x_{\mathrm{out}}t`$ correlation term $`D_{x_{\mathrm{out}},\beta _{}t}`$ in Eq.(5). In the inset of the lower-right panel of Fig. 2, we show the ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ for midrapidity charged pions. It is seen that the ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ obtained from the emission function (solid curve) has a value between 1.0 and 1.3 as in predictions based on the hydrodynamical model with freezeout treated via the hadronic transport model Soff:2000eh . However, with the radius parameters extracted from the Gaussian fit to the correlation function, the ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ becomes much closer to 1, similar to the experimental values extracted from the measured correlation function Adler:2001zd .

To investigate the reason for the large difference between the source radii obtained directly from the emission function and from the Gaussian fit to the correlation function, we show in Fig. 3 the $`x_{\mathrm{out}}x_{\mathrm{side}}`$ distribution (left panel) and the $`x_{\mathrm{out}}t`$ distribution (right panel) at freezeout for midrapidity pions with $`125<p_\mathrm{T}<225`$ MeV$`/c`$ from the AMPT model with string melting and $`\sigma _p=10`$ mb. It is seen that the emission source is shifted in the direction of the pion transverse momentum, i.e, $`x_{\mathrm{out}}>0`$. This positive shift in $`x_{\mathrm{out}}`$ results from the collective expansion of the emission source Lin:2001zk .

If the width in $`x_{\mathrm{out}}`$ were much smaller than the average value $`x_{\mathrm{out}}`$, the emission source would have a shell-like shape, similar to the emission source with strong opacity in the hydrodynamical model Tomasik:1998qt . The emission source also shows a large halo around a central core. The halo consists not only of pions from decays of long-lived resonances such as the $`\omega `$ but also of thermal pions. In calculating the correlation function, we have included pions from the decay of $`\eta `$ resonances as in experiments. Their effects on the radius parameters obtained from the Gaussian fit to the correlation function are thus included. On the other hand, we have excluded these pions in evaluating the source radii from the emission function due to the long lifetime of $`\eta `$. Since long-lived resonances mainly affect the correlation function at small relative momenta Csorgo:1994in , they are important in determining the $`\lambda `$ parameter.

As to the $`x_{\mathrm{out}}t`$ distribution of the emission function shown in the right panel of Fig. 3, it has a strong positive $`x_{\mathrm{out}}t`$ correlation as clearly seen from the solid curve with open diamonds, which shows that the average value $`x_{\mathrm{out}}`$ increases with the freezeout time $`t`$. The $`x_{\mathrm{out}}t`$ correlation leads to a large positive value for the $`x_{\mathrm{out}}t`$ correlation term $`D_{x_{\mathrm{out}},\beta _{}t}`$ in Eq.(5). For pions included in generating Fig. 3, the value of $`D_{x_{\mathrm{out}},\beta _{}t}`$ is 168 fm2 and is appreciable compared to 185 and 431 fm2, respectively, for the first and last terms in Eq.(5). This makes it difficult to extract information about the duration of emission from the ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$. Similar results have been seen previously in studies based on the RQMD model at SPS Sullivan:wb ; Fields:sj . We note that the imaging method imaging , developed for extracting the emission function of a source from the correlation function, will be very useful for verifying the non-Gaussian features of the emission source in high energy heavy ion collisions.

In summary, we have studied the pion interferometry in relativistic heavy ion collisions at RHIC. Using a multiphase transport model that includes both initial partonic and final hadronic interactions, we find that the two-pion correlation function is sensitive to the parton scattering cross section, which controls the density at which the parton-to-hadron transition occurs in the AMPT transport model. To reproduce the measured correlation function in central Au+Au collisions at $`\sqrt{s}=130A`$ GeV requires both the melting of initial strings to partons and a large parton scattering cross section, i.e., a low parton density for the transition from the partonic to the hadronic matter. We further find that the emission source is non-Gaussian in space and time. It not only shifts significantly to the direction along the pion transverse momentum but also has a strong correlation between this displacement and the freezeout time. Consequently, the source radii extracted directly from the emission function are about a factor of 2 to 3 larger than the radius parameters extracted from a Gaussian fit to the three-dimensional correlation function. Furthermore, the ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ obtained from the emission function is larger than that extracted from a Gaussian fit to the correlation function, which is found to be close to one. Although the correlation function requires only the space-time information of pions at freezeout, it is shown to be sensitive to the partonic dynamics during the early stage of heavy ion collisions. The study of pion interferometry thus helps to confirm the formation of the partonic matter at RHIC and to study its properties.

###### Acknowledgements.

We appreciate useful discussions with L.-W. Chen, P. Danielewicz, M. Gyulassy, U. Heinz, M. Lisa, M. Murray, P. Philip, T. Csörgö, S. Pratt, N. Xu, and Q.H. Zhang. This paper is based on work supported by the U.S. National Science Foundation under Grant Nos. PHY-9870038 and PHY-0098805, the Welch Foundation under Grant No. A-1358, and the Texas Advanced Research Program under Grant No. FY99-010366-0081.


text:
# Circumbinary Disk evolution

## 1 Introduction

Observations of main sequence stars have shown that about 60% belong to binary or multiple systems (Duquennoy & Mayor 1991). However, recent observations of pre-main sequence stars seem to indicate that nearly all stars are born in multiple systems (Reviews in Mathieu et al. 2000; Zinnecker & Mathieu 2001).

Apart from many spectroscopic binary systems with surrounding disks like DQ Tau and AK Sco, there have been only two observations of directly imaged wide optical binary systems with circumbinary disks, namely GG Tau (Dutrey et al. 1994) and UY Aur (Duvert et al. 1998). The observations of the spectroscopic binaries usually include emission data from the infrared and optical part of the spectrum. Observations of the luminosities of the boundary layer are typically used to estimate accretion rates onto the binary stars (Gullbring et al. 1998; Hartigan et al. 1995).

Bate & Bonnell (1997) model the early phase of a molecular cloud collapsing onto a protobinary system with an SPH code and derive properties of the created circumbinary and circumstellar disks. We are interested in the following phase of the evolution when the envelope has been cleared. By modeling the evolution we compare fully developed circumbinary and circumstellar disks and their properties with observational data.

In a system consisting of a stellar binary surrounded by an accretion disk there is a continuous exchange of angular momentum between the binary and the disk. As the binary is orbiting with a shorter period than the disk material it is generating positive angular momentum waves in the disk which, upon dissipating, deposit their energy and angular momentum in the disk. Hence, there is a continuous transfer of angular momentum from the binary to the outer disk. Upon receiving this angular momentum, the inner disk speeds up, and the material is receding from the binary, leaving a region of very low density around the central binary, a so called inner gap. At the same time, the loss of angular momentum leads to a secular shrinkage of the semi-major axis of the binary. The width and form of the gap depend on disk parameters such as temperature and viscosity and on the orbital parameters (eccentricity) of the binary as well (Artymowicz et al. 1991; Artymowicz & Lubow 1994).

First calculations dealing with the evolution of circumbinary disks so far have been using the numerical method of smoothed particle hydrodynamics. Applying such a scheme Artymowicz & Lubow (1996) could demonstrate that for sufficiently high viscosities and temperatures material from the disk is able to penetrate into the inner gap region of the disk to become eventually accreted by one of the binary stars. This process which is facilitated by a gas flow along the saddle points of the potential leads to a pulsed accretion whose magnitude and phase dependence for eccentric binaries was estimated. However, when using only a small number of particles all having the same mass only a limited resolution is achievable.

Later Rozyczka & Laughlin (1997) used a finite difference method in attacking this problem. Using a constant kinematic viscosity coefficient and eccentric binaries they also found, in agreement with Artymowicz & Lubow (1996), a pulsed accretion flow across the gap, with most of the accretion occurring onto the lower mass secondary star.

In this paper we extend these computations and solve explicitly a time dependent energy equation which includes the effects of viscous heating and radiative cooling. We aim at the structure and dynamics of the disk as well as the gas flow in the close vicinity of the binary star. To that purpose we utilize a newly developed method which enables us to cover the whole spatial domain. For the first time we perform long-time integration of the complete system covering several hundred orbital periods of the binary and compare the properties of the evolved systems with observational data such as spectral energy distributions in the infrared and optical bands and accretion rates estimated from luminosities.

In the next section we layout the physical model followed by some remarks concerning important numerical issues (Sect. 3). We describe the setup of the various numerical models in Sect. 4. The main results are presented in Sect. 5 and our conclusions are given in Sect. 6.

## 2 Physical Model

### 2.1 General Layout

The system consists of a stellar binary system and a surrounding circumbinary accretion disk. We assume a geometrically thin disk and describe the disk structure by means of a two-dimensional, infinitesimal thin model in the $`z=0`$ plane with the origin at the center of mass (COM) of the binary. We use vertically averaged quantities, such as the surface mass density

$$\mathrm{\Sigma }=_{\mathrm{}}^{\mathrm{}}\rho 𝑑z,$$

where $`\rho `$ is the regular three-dimensional density. We work with two different coordinate systems. Firstly, a cylindrical coordinate system $`(r,\phi )`$ centered on the COM of the binary stars. As such a coordinate system does not allow for a full coverage of the plane, we overlay the center additionally with a Cartesian $`(x,y)`$ grid (see Sec. 3.2).

The gas in the disk is non-self-gravitating and is orbiting a binary system with stellar masses $`M_1`$ and $`M_2`$ which are modeled by the gravitational potential of two solid spheres. The binary is fixed on a Keplerian orbit with given semimajor axis $`a`$ and eccentricity $`e`$. Hence, we neglect self-gravity and the gravitational backreaction of the disk onto the stars. This is justified because the total mass of the disk $`M_\mathrm{d}`$ within the simulated region, which extends up to several tens of semi-major axis $`a`$ of the binary, is typically much smaller than that of the central binary. Hence, the timescale for a change in $`a`$ and $`e`$ is indeed much longer than the timescales considered here (Artymowicz et al. 1991). For the majority of the computations we employ an inertial coordinate system, however test calculations were performed in a system corotating with the binary as well.

### 2.2 Equations

The evolution of the disk is given by the two-dimensional ($`x,y`$ or $`r,\phi `$) evolutionary equations for the density $`\mathrm{\Sigma }`$, the velocity field $`u(u_x,u_y)=(u_r,u_\phi )`$, and the temperature $`T`$. In a coordinate-free representation the equations read

$$\frac{\mathrm{\Sigma }}{t}+\left(\mathrm{\Sigma }u\right)=0$$
(1)
$$\frac{\mathrm{\Sigma }u}{t}+\left(\mathrm{\Sigma }uu\right)=P\mathrm{\Sigma }\mathrm{\Phi }+𝐓$$
(2)
$$\frac{\mathrm{\Sigma }\epsilon _{\mathrm{tot}}}{t}+\left[\left(\mathrm{\Sigma }\epsilon _{\mathrm{tot}}+P\right)u\right]=\mathrm{\Sigma }u\mathrm{\Phi }\left(Fu𝐓\right)$$
(3)
Here $`\epsilon _{\mathrm{tot}}=1/2u^2+\epsilon _{\mathrm{th}}`$ is the sum of the specific internal ($`\epsilon _{\mathrm{th}}=c_\mathrm{v}T`$) and kinetic energy. $`P`$ is the vertically integrated (two-dimensional) pressure

$$P=\frac{R\mathrm{\Sigma }T}{\mu }$$

with the midplane temperature $`T`$ and the mean molecular weight $`\mu `$, which is obtained by solving the Saha rate equations for Hydrogen dissociation and ionization and Helium ionization.

The gravitational potential $`\mathrm{\Phi }`$, generated by the binary stars, is given by

$$\mathrm{\Phi }=\underset{i=1,2}{}\{\begin{array}{cc}\frac{GM_i}{\left|rr_i\right|}\hfill & \mathrm{for}|rr_i|>R_{}\hfill \\ \frac{GM_i\left(3R_{}^2\left|rr_i\right|^2\right)}{2R_{}^3}\hfill & \mathrm{for}|rr_i|R_{}\hfill \end{array}$$
(4)
where $`r_i`$ are the radius vectors to the two stars, and $`R_{}`$ is the stellar radius assumed to be identical for both stars. The second case in Eq. (4) gives the potential inside a star with radius $`R_{}`$ having a homogeneous (constant) density.

### 2.3 Viscosity and Disk Height

The effects of viscosity are contained in the viscous stress tensor $`𝐓`$. Here we assume that the accretion disk may be described as a viscous medium driven by some internal turbulence which we approximate with a Reynolds ansatz for the stress tensor. The components of $`𝐓`$ in different coordinate systems are spelled out explicitly for example in Tassoul (1978). A useful form for disk calculations considering angular momentum conservation is given in Kley (1999).

For the kinematic shear viscosity $`\nu `$ we assume in this work an $`\alpha `$-model (Shakura & Sunyaev 1973) in the form

$$\nu =\alpha c_\mathrm{s}H$$
(5)
where $`c_\mathrm{s}`$ and $`H`$ are the local sound-speed and vertical height, respectively. The local vertical height $`H(r)`$ is computed from the vertical hydrostatic equilibrium

$$\frac{\mathrm{\Phi }}{z}=\frac{1}{\rho }\frac{p}{z}$$
(6)
where $`\rho `$ and $`p`$ are the standard three dimensional density and pressure. Substituting Eq. (4) for the potential with $`|rr_i|>R_{}`$, assuming an ideal gas law $`p=R\rho T/\mu `$, and integrating over $`z`$ yields

$$\rho =\rho _0e^{\left(\frac{1}{2}\frac{z^2}{H^2}\right)}$$
(7)
where $`\rho _0=\mathrm{\Sigma }/[(2\pi )^{1/2}H]`$ is the midplane density, and $`H`$ is the vertical height given by

$$H(r)=\left(\underset{i=1,2}{}\frac{GM_i}{c_\mathrm{s}^2\left|rr_i\right|^3}\right)^{\frac{1}{2}}=\left(\underset{i=1,2}{}H_i^2(r)\right)^{\frac{1}{2}}$$
(8)
which can be split into the single star disk heights given by $`H_i(r)=c_\mathrm{s}\sqrt{\frac{\left|rr_i\right|^3}{GM_i}}`$. The sound speed is given here by $`c_\mathrm{s}=RT/\mu `$. We assume a vanishing physical bulk viscosity $`\zeta `$, but consider it in the artificial viscosity coefficient (see Kley 1999).

### 2.4 Radiative Balance

The influence of radiative and viscous effects on the disk temperature are treated in a local balance. The balance equation for the internal and radiative energy considering only these two effects reads

$$\frac{\left(\mathrm{\Sigma }c_\mathrm{v}T+aT^4\right)}{t}=Q_{\mathrm{diss}}Q_{\mathrm{rad}}$$
(9)
where $`Q_{\mathrm{diss}}`$ and $`Q_{\mathrm{rad}}`$ denote the viscous dissipation and radiative losses, respectively. For the viscous losses we use the vertically averaged expression

$$Q_{\mathrm{diss}}=u𝐓=\frac{1}{2\mathrm{\Sigma }\nu }Tr\left(𝐓^2\right)$$
(10)
For the radiative transport

$$Q_{\mathrm{rad}}=F_0\frac{F_z}{z}$$
(11)
where $`F_0`$ is the flux vector in the $`z=0`$ plane. Here we consider only the losses in the vertical direction, i.e. $`F_0=0`$, a standard approximation in accretion disk theory. Integration over the vertical direction yields

$$Q_{\mathrm{rad}}=2F_{\mathrm{rad}}=2\sigma _\mathrm{B}T_{\mathrm{eff}}^4$$
(12)
with the local effective (surface) temperature $`T_{\mathrm{eff}}`$ which is related to the midplane temperature $`T`$ through

$$T_{\mathrm{eff}}=\tau ^{1/4}T=\left(\frac{1}{2}+\frac{3}{4}\kappa \mathrm{\Sigma }\right)^{1/4}T$$
(13)
using an interpolated opacity $`\kappa (\rho _0,T)`$ adapted from Lin & Papaloizou (1985). Note that this formulation of radiative transfer is a very good approximation only in the optically thick regime.

## 3 Numerical Issues

In order to study the disk evolution, we utilize a finite difference method to solve the hydrodynamic equations outlined in the previous section. As we intend to model possible accretion flow onto the binary stars and to achieve a very good resolution in their vicinity, we use a Dual-Grid-technique.

### 3.1 General

The code is based on the hydrodynamics code RH2D suited to study general two-dimensional systems including radiative transport (Kley 1989). RH2D uses a spatially second-order accurate, mixed explicit and implicit method. Due to an operator-splitting method it is semi-second order accurate in time. The advection is computed by means of the second order monotonic transport algorithm, introduced by van Leer (1977), which guarantees global conservation of mass and angular momentum. Advection and forces are solved explicitly, and the viscosity and radiation are treated implicitly. The formalism for application to thin disks in $`(r\phi )`$ geometry has been described in detail in (Kley 1998, 1999). Here this scheme is extended to allow for a Dual-Grid system (see 3.2) and a more detailed energy balance Eq. (9).

As the calculations have shown, the system behaves extremely dynamically with very strong, moving gradients at the inner edge of the disk. To model such a situation, the viscosity has to be treated implicitly to ensure numerical stability. The coupling of $`u_r`$ and $`u_\phi `$ requires a solution method similar to that in Kley (1989). For stability purposes an additional implicit artificial viscosity has to be included in the computations (Kley 1999).

This code has been employed successfully in related simulations dealing with embedded planets in disks (Kley 1999, 2000). The numerical details of the present simulations are presented in Günther (2001).

### 3.2 Dual-Grid technique

The Dual-Grid technique employed here combines two desirable features of thin disk computations: a) An ideally suited coordinate system $`(r\phi )`$ in the outer bulk parts of the disk which ensures conservation of angular momentum, and b) Coverage of the region near the central binary star allowing the determination of the mass flow onto the stars.

In the main part the domain is covered by a cylindrical grid. On this (outer) grid the equations are formulated in polar coordinates, and the solution technique guarantees exact conservation of global angular momentum. The conservation property is absolutely crucial in obtaining physically reliable results for long term disk evolution (Kley 1998).

However, such a grid cannot be extended to the very center and leaves a hole in the middle. For typical computations this does not pose serious problems, but in this case we are dealing with eccentric binary systems in the middle. Hence, the stars or their circumstellar disks may easily cross the inner grid boundary causing strong irregularities. In any case a central hole does not allow for an accurate determination of the individual mass accretion rates onto or a possible overflow between the stars. Therefore, using a Cartesian grid centered on the origin enables an accurate calculation of the mass flows in the vicinity of the stars. Applying this technique, a high local resolution can be obtained in the center as well. An example of such a grid structure is displayed in Fig. 1.

Related numerical schemes using a nested grid technique, though not in different coordinate systems, have been applied in astrophysical simulations by a number of authors (eg. Ruffert 1992; Yorke et al. 1993).

The necessary equations (in Cartesian and cylindrical coordinates) are then integrated, independently and in parallel, on the two grids, and after each timestep the necessary information in the overlap region has to be exchanged. Restrictions are imposed on the time step for stability reasons, the Courant-Friedrichs-Lewy (CFL) condition must be fulfilled during each integration, on each grid.

Test calculations using a single Cartesian grid covering the whole domain have shown that the non-conservation of angular momentum yields serious deviations from known analytical results in particular for longer integration times.

As the grids possess a different geometry, exact conservation of momenta and energy cannot be guaranteed in the conversion/interpolation process which uses linear interpolation to exchange density, temperature and velocity. Exchanging density, internal energy and momentum was also implemented and leads to more correct conservation of momenta and energy. However, the advantage of a complete coverage of the physical space in between the two stars clearly outweighs this disadvantage.

### 3.3 Tests

We ran several tests to check the accuracy and numerical stability of the grid overlap region, including an infalling shock and a Sedov like explosion. For these tests no viscosity nor radiative effects are applied. All tests show that the grid overlap technique is sufficient for handling problems with radial symmetric features such as accretion disks.

The shock experiment initial setup is a step function of density and temperature with the floor starting at $`r_0`$ (about 10 grid cells outside of the inner boundary of the cylindrical grid) and covering the whole Cartesian grid:

$$\mathrm{\Sigma }(r)=\mathrm{\Sigma }_{\mathrm{min}}+\mathrm{\Sigma }_0\mathrm{\Theta }(rr_0).$$
(14)
In Fig. 2 you can see the grid setup and the shock in the initial infalling phase just after crossing the grid boundaries. The shock front remains radially symmetric and no reflections are visible. After the shock compresses into a single point it switches to the outflowing phase. You can see the shock front propagating outwards just after crossing the grid boundaries in Fig. 3. Apart from unavoidable grid effects near the center the shock front remains nicely circular and again neither disturbances nor reflections are visible at the grid interface.

## 4 General model design

The main goal of this study is the investigation of the characteristic features of the hydrodynamic flow of the circumbinary disk.

From now on, we refer to non-dimensional units. All the lengths are expressed in units of the semi-major axis $`a`$ of the binary. This is constant because the time scale on which the binary orbital elements change due to mass and angular momentum transfer is much longer than the dynamical time scales considered here (see Artymowicz & Lubow 2001). Masses are in units of the solar mass and time is given in units of the binary’s orbital period.

The whole azimuthal range of the disk is taken into account by considering a computational domain represented by $`2\pi \times [0,r_{\mathrm{out}}]`$, where typically $`r_{\mathrm{out}}=10`$ in units of the binary separation. The cylindrical grid covers a region $`2\pi \times [r_{\mathrm{in}},r_{\mathrm{out}}]`$ with $`r_{\mathrm{in}}=0.1`$ represented by up to $`256\times 256`$ meshpoints, and the central Cartesian grid has up to $`128\times 128`$ meshpoints covering $`2r_{\mathrm{in}}\times 2r_{\mathrm{in}}`$ with an additional small overlap.

### 4.1 Accretion onto the Stars

As will be seen, matter flows from the disk through the gap and circulates the stars in small circumstellar disks. From those it can be accreted onto the stars. This stellar accretion is accounted for by removing some mass from a region close to the stars. This removal is performed on both grids dependent on the star positions. The removed mass is not added to the dynamical mass of the stars, but is just monitored.

The region $`R`$ out of which matter is removed is given as a fraction (0.2) of the Roche lobe size for the wide systems and as the stellar cores for the tight spectroscopic systems. The accretion process typically involves only a few grid cells on each side of the stars for the wide systems and the whole resolved stellar core (up to 11x11 grid-cells) for the tight systems, making it a locally confined process.

Two methods for deciding whether a grid cell is targeted for accreting a part of its mass were implemented. Accreting a constant fraction $`f`$ per time of the mass exceeding a defined minimum density $`\mathrm{\Sigma }_{\mathrm{min}}`$ leads to a very smooth accretion process. The accreted mass per iteration is computed by

$$M_{\mathrm{acc}}=\underset{i,jR}{}\mathrm{min}(1.0,f\mathrm{\Delta }t)\mathrm{max}(0.0,\mathrm{\Sigma }_{ij}\mathrm{\Sigma }_{\mathrm{min}})\mathrm{\Delta }\mathrm{Vol}_{ij}$$
(15)
Limiting the density in the accretion region $`R`$ to an arbitrary value $`\mathrm{\Sigma }_{\mathrm{max}}`$ results in a more realistic structure of the circumstellar disks and to more correct spectral energy distributions of the circumstellar disks.

$$M_{\mathrm{acc}}=\underset{i,jR}{}\mathrm{max}(0.0,\mathrm{\Sigma }_{ij}\mathrm{\Sigma }_{\mathrm{max}})\mathrm{\Delta }\mathrm{Vol}_{ij}$$
(16)
This accretion process is not smooth, but consists of accretion events because exceeding the maximum density is a discrete process and not distributed evenly over time. Nevertheless the time-averaged accretion rates are the same for both methods.

The inferred accretion rates in any event compare favorably with the average disk accretion rate, as estimated from the radial mass fluxes through different surfaces, and the evolution of the total mass found in the circumstellar disks.

### 4.2 Spectra generation

Spectral energy distributions are generated decoupled from the dynamic evolution of the disk as a post-processing step. Here the model from Adams et al. (1988) is extended to work on data generated by our numerical simulation. The flux at the observer is obtained by integrating over the disk surface assuming local black-body radiation

$$F_\nu =\frac{\mathrm{cos}i}{D^2}_{r_0}^{R_d}_0^{2\pi }B_\nu \left[T(r,\phi )\right]\left(1e^{\tau (r,\phi )}\right)rd\phi dr$$
(17)
where $`i`$ is the inclination of the disk, $`D`$ the distance to the observer and $`\tau `$ the line-of-sight optical depth through the disk. $`\tau `$ can be obtained using the opacity $`\kappa `$, $`\tau (r,\phi )=\frac{\kappa \mathrm{\Sigma }(r,\phi )}{\mathrm{cos}i}`$. $`T`$ is the surface disk temperature which has to be estimated using the optical depth of the disk obtained from an opacity table.

The integration is limited to regions outside of the stellar cores as temperatures and densities from inside the cores can be nothing but wrong. The star emission is accounted for by adding a black-body spectrum from an appropriate effective star temperature.

### 4.3 Initial and boundary conditions

The initial density distribution is proportional to $`r^d`$, where $`d`$ is either set to $`0.5`$ which is the equilibrium density distribution for an accretion disk around a central star having a constant kinematic viscosity, or given by a fit of the generated spectral energy distribution to the one observed. However, we superimpose to this an axisymmetric gap around the COM of the two stars, given by an approximate gap radius $`r_{\mathrm{gap}}`$ (Artymowicz & Lubow 1994) and a gap function

$$f_{\mathrm{gap}}:=\left(e^{\frac{rr_{\mathrm{gap}}}{0.1r_{\mathrm{gap}}}}+1\right)^1.$$
(18)
Figure 4 shows a typical surface density at $`t=0`$. The initial velocity field is that of a Keplerian disk orbiting the center of mass of the two stars.

The gravitational potential in Eq. (4) is phased-in smoothly from an initial central potential for the total mass, $`M_1+M_2`$, during an initialization phase lasting typically ten orbital periods. This allows the system to adjust smoothly from the semi-stationary setup of a single central star to the binary star situation.

For the boundary conditions, periodicity is imposed at $`\phi =0`$ and $`\phi =2\pi `$. We set outflowing boundary conditions at the outer radial border ($`r_{\mathrm{out}}`$), i.e. the radial velocity gradient at $`R_{\mathrm{max}}`$ is zero, and truncated to positive values. Angular velocity there equals the unperturbed initial Keplerian value.

For single grid calculations with a cylindrical grid we use zero gradient boundary conditions for both velocity components and the density at the inner border.

## 5 Results

Before discussing the results concerning particular astronomical objects, we would like to present the details of a specific test case of an equal mass binary star on a circular orbit. Such a system is well suited to check important properties (symmetries, stability) of the numerical method used. Also general properties of accretion disks in binary systems can be derived from such a test case.

### 5.1 Equal Mass Binary

The equal mass binary on a circular orbit serves as a testing platform for numerical stability and accuracy. In a coordinate frame corotating with the binary one expects after some transient time a quasi-stationary situation, which is symmetric with respect to the line connecting the stars. As the viscosity is solved implicitly by solving a large matrix system iteratively, exact symmetry cannot be preserved numerically. Also usual density contrasts are of the order of $`10^3`$ for numerical floor to circumstellar disks and from circumstellar disks to the circumbinary disk. Hence, this test serves as an excellent example for estimating the ability to preserve symmetry over long time scales.

In Fig. 5 you can see a circumbinary disk after 90 orbits evolution time in almost perfect quasi-stationary state. The deviation from line-symmetry is small. Only with very high resolution simulations we observe spiral features in the outer regions of the disk, while spiral arms originating from the inner gap of the circumbinary disk flowing onto the circumstellar disks can be observed even for low resolutions. Those spiral arms feeding the circumstellar disks remain stable for the whole period of a non-eccentric binary.

The accretion rates for both accretion mechanisms onto the two stars are nearly identical and do not show any significant variations for different numerical parameters. Thus, our implementation of the accretion mechanism is robust.

Changing model parameters such as $`\mathrm{\Sigma }_0`$, gap width, or disk size we obtain characteristic variations of the resulting spectral energy distribution of an evolved disk. Figures 6-8 show dependencies of these parameters on the spectral energy distribution that are used to fit parameters of the spectroscopic systems to match observed spectra.

### 5.2 Spectroscopic Systems

Hereafter we discuss models for the two spectroscopic young binary stars DQ Tau and AK Sco which have a separation of only a fraction of an $`\mathrm{AU}`$. The physical parameter for the systems are taken from Andersen et al. (1989); Favata et al. (1998); Jensen & Mathieu (1997) for AK Sco and Mathieu et al. (1997) for DQ Tau. The relevant information for the simulations has been summarized in Table 1, disk masses $`M_\mathrm{d}`$ are integrated densities over the computational domain specified by the disk radii $`R_\mathrm{d}`$. Images showing the disk configuration in apastron and periastron phase are shown in Figs. 9 and 10.

Generally, for eccentric systems we observe spiral features in the outer circumbinary disk only for very high resolution simulations, while again spiral arms feed the circumstellar disks even with low resolutions. Opposed to the non-eccentric systems the spiral arms develop during apastron phase and are torn off together with the outer parts of the circumstellar disks after periastron phase for high eccentric systems. This tearing off is caused by raising centrifugal forces acting on the circumstellar material.

For both systems we compared the accretion rates resulting from the simulation with accretion rates derived from of observations as done by Hartigan et al. (1995) and Gullbring et al. (1998) (see Table 2). As no data was available for AK Sco, a comparison with observational data was not possible, but as both spectroscopic systems have similar system parameters, accretion rates of the same order of magnitude can be expected.

As can be seen from Figs. 11 and 12 the accretion process happens synchronized with periastron phase of the binaries. This is in agreement with numerical simulations from Artymowicz & Lubow (1996) and Rozyczka & Laughlin (1997) and with observations which show excess luminosity at these phases (Mathieu et al. 1997).

From our simulations we can derive an averaged accretion rate of $`0.8310^8\mathrm{M}_{}\mathrm{yr}^1`$ for the primary, $`0.8210^8\mathrm{M}_{}\mathrm{yr}^1`$ for the secondary of AK Sco and $`0.6210^8\mathrm{M}_{}\mathrm{yr}^1`$ for the primary, $`0.5110^8\mathrm{M}_{}\mathrm{yr}^1`$ for the secondary of DQ Tau. These accretion rates are in good agreement with the observational rates.

It is expected that the luminosity increases in periastron because tidal compression of the disk leads to a temperature rise with subsequent luminosity increase. We tried to generate luminosity curves to show that indeed excess luminosity is produced at periastron phase, but unfortunately those curves show both brightening and darkening dependent on the implemented accretion mechanism. At the very low averaged accretion rates of $`\dot{M}710^9\mathrm{M}_{}\mathrm{yr}^1`$ the disk is very optically thin. In that regime the implemented radiative loss mechanism is not very accurate and thus these results are not meaningful.

For the spectroscopic binaries we can generate spectral energy distributions out of the simulation data which can be made to match the observed ones by manually adjusting initial model parameters such as disk mass, gap width and density profile $`d`$. For DQ Tau Fig. 13 shows the spectral energy distribution composed out of the disk and the star, Fig. 14 shows the one for the AK Sco system. The effective stellar temperatures used in the model have been adjusted to match the observational data. They differ from the quoted values in the literature in Table 1. This may be due to difficulties in separating individual stellar and disk contributions to the total flux in models and observations. Also occultation is not accounted for on adding the star spectrum to the disk.

### 5.3 Wider Systems

There are only two well known systems where the binary and the disk have been imaged directly, GG Tau (Dutrey et al. 1994; Guilloteau et al. 1999) and UY Aur (Close et al. 1998; Duvert et al. 1998). These systems have a much wider separation of tens to over 100 AU, a mass ratio different from unity and show consequently a quite different behavior. In Table 3 the system parameters used for the simulations are shown.

From our simulations we derive averaged accretion rates for both UY Aur and GG Tau that are comparable with the observed quantities in Table 2. In the case of UY Aur we get $`0.510^7\mathrm{M}_{}\mathrm{yr}^1`$ for the primary and $`0.2510^6\mathrm{M}_{}\mathrm{yr}^1`$ for the secondary, respectively. Simulations of GG Tau lead to accretion rates of $`0.1110^8\mathrm{M}_{}\mathrm{yr}^1`$ for the primary and $`0.3310^9\mathrm{M}_{}\mathrm{yr}^1`$ for the secondary. As expected from the difference of the disk masses ($`0.15\mathrm{M}_{}`$ vs. $`1.1\mathrm{M}_{}`$) the accretion rate for GG Tau is one order of magnitude less than the one for UY Aur.

Accretion is reinforced at periastron phases as for the spectroscopic systems, but events are an order of magnitude lower for GG Tau and nearly vanish for UY Aur due to the lower eccentricity of the systems. Also both systems show more phase dependent accretion for the secondary than for the primary.

In Figs. 15 and 16 you can see the inner parts of the simulated circumbinary disks. The circumstellar disks remain intact over the whole period. It was not possible to generate meaningful spectral energy distributions of the disks as they are too cold for these wide systems. Also observational data is lacking for these systems.

## 6 Conclusion

With the present calculations we have modeled the accretion flow onto the stars within a circumbinary disk in more detail using more realistic physics, much higher resolution and a much longer time integration. Our main achievements may be summarized as follows:

1. Development of a new Dual-Grid technique combining two coordinate systems, which is ideally suited for modeling planar accretion disks including their central objects.
2. Performing long-time integrations of binary stars and circumbinary disks, including a detailed vertical energy balance.
3. For close binaries with separations of only a fraction of an AU we find a relatively large accretion rate crossing the gap of the order of $`10^8\mathrm{M}_{}\mathrm{yr}^1`$.
4. For wide binaries $`(a=10200\mathrm{AU})`$ the accretion rate is comparable to those of the close binaries only because of the large disk masses of the two systems.
5. For a suitable variation of the initial density distribution $`\mathrm{\Sigma }r^d`$ of the disk, with $`d`$ between $`1.5`$ and $`2`$ the spectral energy distributions of DQ Tau and AK Sco could be fitted well to the observational data.

Future models of this kind need to include a more detailed structure of the stars as well as a more detailed model of the accretion process. The numerical resolution we have reached by now is so fine that the individual stars in close binaries are already resolved by usually about 100 grid-cells within low-res calculations.

To obtain a detailed comparison with observed light curves in particular the phase dependence a more elaborate model of the radiative loss of the very thin material in the circumstellar disks is required. In that case 3d simulations may be necessary which are beyond present day computational resources.

###### Acknowledgements.

Part of this work was funded by the German Science Foundation (DFG) under SFB 382 Simulation physikalischer Prozesse auf Höchstleistungsrechnern. We would like to thank Dr. Hubert Klahr for many stimulating discussions during this project.