Femtoscopy in Relativistic Heavy Ion Collisions and its Relation to Bulk Properties of QCD Matter
Abstract
Using a viscous hydrodynamic model coupled to a hadronic cascade code, numerous features of the dynamics and equilibrium properties are explored for their impact on femtoscopic measurements. The equation of state, viscous parameters and initial conditions are investigated. We find that femtoscopy is affected by numerous model features at the 10% level, and that by including features and adjusting unknown parameters, one can explain experimental source size measurements to better than 10%.
pacs:
25.75.Gz,25.75.LdI Overview
The bulk properties of QCD matter, as created in relativistic heavy ion collisions, largely manifest themselves in soft hadronic observables of particles with transverse momentum less than one GeV/. These observables can be divided into three classes: spectra, flow (or largescale correlations), and correlations at small relative momentum. This last class is referred to as femtoscopy Lisa:2005dd () since these correlations are used to determine spacetime characteristics of emitting sources. Correlation functions, , can be linked to the outgoing phasespace distributions, or more precisely the source function , which describes the probability that two particles with the same velocity, whose total momentum is , are separated by in their asymptotic trajectory. Due to their inherent sixdimensional nature, correlations have proven to be the most difficult of all RHIC observables to fit with full dynamic models. The measurements are amenable to being fit by simple geometric models of the final state, provided that the models incorporate strong radial collective flow, and a rapid dissolution into a thermal assortment of resonances Kisiel:2008ws (); Kisiel:2006is (); Retiere:2003kf (); Helgesson:1997zz (). However, many dynamic models, especially hybrid hydrodynamic/cascade descriptions, lead to more extended emission durations which lead to signicantly different shapes for the outgoing phase space distributions.
The information in correlations are often reduced to Gaussian source parameters, , and , which are functions of the transverse momentum , and describe the shape of the outgoing phasespace distribution of zerorapidity particles with a specific and a specific azimuthal angle. Here, refers to the longitudinal dimension, along the beam axis, describes the outward dimension, parallel to the momentum, and refers to the sideward dimension, perpendicular to the beam and to the particle’s velocity. Asymptotically, the Gaussian form fits the phase space density to the form,
(1) 
The term is irrelevant for identical particles since correlations are only sensitive to the relative position of two particles of the same velocity. For nonidentical particles, one would also be sensitive to the relative position of centroids for the two species, .
This study focuses on how these parameters are affected by choices made in modeling the reaction. The femtoscopic relation to the equation of state has long been studied. First, a stiffer equation of state leads to more rapid expansions, with emission at earlier times and more confined to a brief burst. The reduction in the mean emission time reduces as the system has less time to expand longitudinally before emission. The increased suddenness leads to a shorter relative to , as longlived emission allows those particles emitted earlier to move ahead of the lateremitted particles along the outward direction. Furthermore, a softer equation of state leads to higher entropy, and for fixed spectra, the entropy will be grow for increasing source volumes, . Since the total entropy can be ascertained in a quasimodelindependent fashion from spectra and source dimensions, and since entropy is conserved during the expansion, the product of the three dimensions is strongly linked to the equation of state Pal:2003rz ().
Femtoscopic source sizes are also affected by nonequilibrium aspects of the dynamics. Bulk viscosity, which is expected to be significant in the neighborhood of the critical temperature when the system struggles to maintain equilibrium, lowers the effective pressure and thus increases the entropy and leads to larger source dimensions. Shear viscosity is mainly important at early times, when velocity gradients are large and highly anisotropic. This leads to an enhancement in the transverse components of the energy tensor, which accelerates the transverse expansion and gives smaller values of and relative to romatschke (); Pratt:2006ss (). At the earliest times, before even viscous hydrodynamics is applicable, the preequilibrium state might be dominated by longitudinal color fields. These fields can, in principle, lead to exceptionally strong transverse components to the stressenergy tensor, which would amplify the effects of shear viscosity. The importance of early acceleration in explaining the experimental ratio has also been studied by incorporating initial transverse flow into ideal hydrodynamics Gyulassy:2007zz () and by adding a strong repulsive potential into microscopic cascade codes Li:2007yd (). Freestreaming during the first fm/ increases the transverse pressure relative to the longitudinal pressure, which increases radial flow more than elliptic flow Heinz:2002rs (). This had been thought to make it difficult to simultaneously fit both spectra and elliptic flow, though this was accomplished in Broniowski:2008vp (). This issue might be resolved by better understanding the interface between the initialstate description and the following hydrodynamic description.
The initial density profile also affects acceleration at early times Broniowski:2008vp (). Within the typical nuclear crosssection (40 mb), a single nucleon will interact with multiple nucleons from the opposite beam. Depending on the theoretical picture, e.g., colorglasscondensateinspired or wounded nucleon, the average radius of the initial fireball can vary by , with a more compact source being more explosive and leading to smaller values of .
To investigate these effects we apply a relativistic viscous hydrodynamic model, which couples to a cascade code for modeling the hadronic breakup stage and the decay of outgoing resonances. The outgoing phase space points for pions are used to generate source functions, and after convoluting with the squared wave function, generate twopion correlations. These are then treated as data and fit to correlation functions from Gaussian sources to extract , and as a function of . The hydrodynamic code uses IsraelStewart equations which are modified to allow one to tune to the anisotropies of the initial stressenergy tensor. Both the hydrodynamic and cascade descriptions are built on an assumption of azimuthal symmetry and boostinvariance. This prohibits a simultaneous analysis of elliptic flow, or a study of longitudinal acceleration, which is known to affect results at the 510% level. In addition to investigating all the sensitivities alluded to above, we compare data from STAR collaboration at the Relativistic Heavy Ion Collider (RHIC). We do find solutions which come close to the the data without employing any particularly disquieting assumptions or any parameters outside what we would consider a reasonable range. Although this paper focuses on femtoscopy, the mean of various calculations is also presented and compared to data.
After reviewing details of the model in the next section, the following section is devoted to the effects of varying the equation of state, viscosities, and initial conditions. The summary is devoted to drawing conclusions with an emphasis on understanding what future improvements in models and analysis are needed to reach rigorous quantitative statements about the microscopic properties of the QCD matter formed in relativistic heavy ion collisions.
Ii The Model
For generating interferometric source functions, phase space points are first generated from a viscous hydrodynamic model, then fed into a cascade model which models the lowdensity hadronic stage of the collision. Both models are written in terms of the variables , , and , where is the proper time, represents the longitudinal position, and and represent the radial position and azimuthal angle. Both models were developed assuming radial symmetry and boost invariance which eliminates and from consideration. By reducing the dimensionality, both speed and accuracy are vastly improved. The viscous hydrodynamic model is based on the formalism in Pratt:2007gj (), with more details provided below. The subsequent subsections provide details of the three model components.
ii.1 Viscous Hydrodynamic Model
First, a review of the modified IsraelStewart formalism described in Pratt:2007gj () is presented. A basic description of viscosity and the NavierStokes equation can be found in weinberg (). Recently, IsraelStewart hydrodynamics israelstewart () has been extended and applied to nuclear physics muronga (); koide (); baierromatschke (); heinzchaudhuri (); songheinz (). In ideal hydrodynamics, the stress energy tensor becomes when viewed in the fluid rest frame (Here latin indices refer to spatial components only). Viscous hydrodynamics deals with the deviation of from . For all the hydrodynamic calculations here, the fluid rest frame is defined such that , and diffusion of conserved particle numbers through fluid elements is ignored. In the fluid frame the deviations of can be expressed through five independent traceless components, , and the deviation of the trace, .
(2)  
The shear components are related on a onetoone basis to the five velocity gradients, ,
(3)  
With these definitions, the NavierStokes equations become
(4) 
For IsraelStewart equations of motion, and are not fixed as is the case for NavierStokes equations, but instead are dynamic objects. The ratios and should decay exponentially toward the NavierStokes values, where and are related to the fluctuation of the stress energy tensor at fixed energy density,
(5)  
where no sum is implied in the last expression. To restrict the values of and to ranges and respectively, the IsraelStewart equations are modified by mapping and to and through hyperbolic tangents. The variables and will follow the IsraelStewart equations above and can become arbitrarily large, while and will be restricted.
(6)  
As derived in Pratt:2007gj () and jou (), the lifetimes, fluctuations and viscosities are not independent,
(7) 
The equations of motion are solved by storing the velocities and energy densities in a mesh defined by the radial coordinate. Following the ideas of the onedimensional calculation in Pratt:2008jj (), mesh points are not stored at equal times but at varying times that enforce local simultaneity, i.e., , where is the four vector describing the separation of two neighboring mesh points. Using the integrated distance along the mesh as seen by comovers,
(8) 
the acceleration in the fluid frame, which equals the rate of change of the transverse rapidity, takes on a simple form,
(9) 
In order to maintain simultaneity between neighboring mesh points, the time step depends on ,
(10) 
To complete the equations of motion, an expression is needed for the evolution of the energy density. This is done by considering the change in the internal energy within a cell defined by adjacent mesh points and a fixed small rapidity range . In the fluid frame, the volume of the cell is
(11) 
where is the radius as viewed in the laboratory frame. In a time step the volume increases both due the increase in the longitudinal dimension and by the increase in the transverse dimensions. Writing
(12)  
the change in the internal energy of the cell is
(13) 
where is the longitudinal direction and refers to the radial direction. Given the internal energy and volume, one then knows the local energy density which closes the equations of motion. For the ideal case, , one recovers , which implies entropy conservation. Indeed, when the code was run in this limit, entropy was conserved to better than 0.2 percent.
The equation of state used for the runs shown here consisted of three parts. For temperatures below 170 MeV, the equation of state was that of a hadronic gas. For a given cell, the pressure was calculated as a function of the energy density and the density of five conserved charges. The conserved charges were the number of strange plus antistrange quarks, the number of baryons plus antibaryons, the effective number of pions (e.g., a meson counts as two pions), the number of s and the number of s. Only the standard octet mesons and octet and decuplet baryons were considered, i.e., the mesons and the baryons. The details of which particle numbers were fixed was not particularly important because the breakup density was chosen to be 400 MeV/fm, which allowed the cells very little time to adjust their chemistry before the evolution was taken over by the cascade code.
For an intermediate range of energy densities, , the equation of state was chosen to have a constant speed of sound, i.e., . Here, is the energy density of an equilibrated hadron gas with a temperature of 170 MeV. In the limit of , the equation of state becomes that of a first order phase transition with latent heat . For energy densities above , the speed of sound was bumped up to to be consistent with lattice gauge theory Cheng:2007jq (). The simple form for the equation of state was used so that by varying and one could study the sensitivity to the equation of state.
The ratio of the shear viscosity to entropy was fixed above . According to the KSS conjecture, this ratio should stay above kss (). Results for varying are investigated in Sec. III. The fluctuation was calculated by considering fluctuations of a massless gas. This gives . The relaxation time from Eq. (7) is then . For , the relaxation time was chosen as , with mb and being the hadron density. This energyaveraged crosssection for a hadronic gas would be somewhat higher than 25 mb, but much of that cross section would be more forward peaked which reduces the effectiveness of collisions to thermalize the matter. It would not be surprising if more sophisticated calculations of the relaxation time would differ by a few tens of percent. The fluctuation was determined by considering the fluctuations inside a hadron gas,
(14) 
The relaxation time is then given by Eq. (7). For the intermediate region, , both and were chosen to vary linearly with the energy density from the hadronic value at to the value for the lower end of the plasma region. Relaxation times were then chosen according to Eq. (7).
Bulk viscosities are expected to be negligible except near Paech:2006st (). This can be understood by considering an isotropically expanding thermalized gas, i.e., a Hubble expansion, which for photons maintains a thermal form to the photon spectrum with the temperature falling . For a nonrelativistic gas, a thermal form also ensues, but with the temperature falling as . If the shape of the momentumspace distribution is already thermal, collisions are needed to maintain the thermal value for . Bulk viscosities thus disappear high above where the temperatures far exceed the quark mass, and for temperatures well below the pion mass.
In contrast, near the degrees of freedom and the condensed fields need to change in order to maintain equilibrium, which leads to a bulk viscosity in that region Paech:2006st (); Karsch:2007jc (). The bulk viscosity was chosen to be zero for and for . In the middle of the mixed region, , was set to a maximum value, . To make the bulk viscosity a continuous function, it was chosen to linearly fall with energy density above and below the maximum value so that it returned to zero at the boundary of the mixed region. Arbitrarily, the relaxation time was chosen to be , which equals the minimum relaxation time for the shear in the plasma phase if one is at the KSS limit. The treatment of the effects of bulk viscosity here are undoubtedly naive. Since the principal source of bulk viscosity might be the nonequilibrium chiral condensate, or field, relaxation times might be very large, the response might be very nonlinear, and the behavior might be oscillatory, which contradicts the IsraelStewart assumption that nonequilibrium deviations decay exponentially. Thus, the investigations can probably only qualitatively describe the impact of nonequilibrium effects towards the trace of the stressenergy tensor. Bulk effects due to nonequilibrium fields would be better treated by a simultaneous solution of the respective wave equations coupled to the hydrodynamic medium.
Figures 1 and 2 display the hydrodynamic evolution for the default parameter set. The kinks in the collective transverse rapidity shown in Fig. 1 arise from the region near . As the matter expands into a region with lower speed of sound a pulse builds up similar to a tsunami. Bulk viscosity, which lowers the effective pressure, , in this region amplifies the pulse. The pulse largely dissipates by the time final breakup occurs as the rapidity profile becomes linear and the energy density profile also becomes smooth.
Viscous effects on the stress energy tensor are illustrated in the three panels of Fig. 2. The effective pressure is displayed in the upper panel. Since the bulk viscosity is set to zero outside the intermediate region, the ratio varies from unity only in this region. If not for the saturation enforced by Eq. (6), the effective pressure might fall below zero. The size of the effect is enhanced by the pulse which results in large velocity gradients at the boundaries of the pulse. The anisotropy of is shown in the lower two panels. At the anisotropy was inserted as a boundary condition with set to zero, which equivalently gives as shown in the middle panel. This is also the saturated value as enforced by Eq. (6), and is maintained for some time due to the large velocity gradients at early times. At the edge of the fireball, where the matter is in the lower density hadronic phase, the strong anisotropy remains due to the large viscosity at low density. However, the behavior in this region is somewhat irrelevant as it is below the breakup density and is handled by the cascade description described in the next section. The lower panel shows , which differs from zero mainly near due to the radial pulse. The large variations shown in Fig. 2 shows how NavierStokes treatments, which can lead to arbitrarily large deviations of the stressenergy tensor, are questionable at early times or in the region of the radial pulse.
The hydrodynamic module was run until the entire system fell below the breakup density. A sampling of emitted hadrons was generated from the entire evolution. At each time step, particles were generated from thermal surface emission of the outermost cell whose energy density was above the breakup density. The generation of particles was consistent with the temperature, density, and the anisotropy of the stressenergy tensor. When a cell’s energy density fell below the breakup density, particles were emitted from that cell in the same manner, except according to volume emission. In order to make such emission consistent with the anisotropy of the stress energy tensor, the particles were first generated according to an isotropic thermal distribution. The momenta , and as seen in the fluid frame were then scaled by factors , and respectively, where
(15) 
The same mechanism is used for surface emission, along with an additional factor taking into account the rate at which the particles leave the surface. For a given particle moving along a collisionless trajectory, the momenta as measured in the local rest frame should scale inversely with the time between collisions. For each component of the momentum, , where is the inverse velocity gradient at the time of the last collision along the given direction and is the time since the last collision. For nonrelativistic particles one can derive the simple scaling form for the various momenta shown above. However, for lighter particles the simple scaling is only approximately justified. Future versions of the program will apply a more sophisticated mechanism for generating particles.
ii.2 Hadronic cascade Model
For energy densities below 400 MeV/fm, a cascade code is used to describe the evolution. The cascade simulates the evolution of the particles as straight line trajectories, punctuated by collisions whose probability is determined by a combination of a fixed cross section of 15 mb, along with resonant absorptions and decays. The resonant cross sections use a simple BreitWigner form with fixed lifetimes, and all collisions and decays are treated as waves. Only resonances from the standard meson octets or from the baryon octet and decuplet are included. This is a simple treatment, with no mean field or Bose effects, but should provide a sufficiently reasonable description of the breakup stage for the interferometric studies presented here.
Particles are entered into the cascade description from a Monte Carlo list generated by the hydrodynamic module. Along with the list of particles, the cascade module is also given a description of the position of the emitting surface as a function of time. Any particle that returns to the interior of the surface during the cascade description is deleted. Assuming that the hydrodynamic code, with its inclusion of shear anisotropies, accurately models the behavior of a hadron gas for energies near the 400 MeV/fm, this should provide a fully consistent interface. For all of the parameter sets studied here, the hypersurface from which particles are created by the hydrodynamic module rapidly collapses resulting in timelike emission for the vast majority of particles. The percentage of particles that are reabsorbed into the hypersurface during the cascade is only about one percent.
One potential issue with many cascade codes is that the finite interaction range, i.e., particles collide at a finite interaction range of approximately , leads to viscous effects Cheng:2001dz (). Usually, such effects are minimized by oversampling the distribution by a factor which reduces the cross sections by and the interaction ranges by . However, in the cascade description applied here, the interaction range is set to zero, thus eliminating such numerical viscosities. This is accomplished by exploiting boost invariance and azimuthal invariance, which allows the trajectories to be treated as radii evolving as a function of the proper time, . When two particles have the same radius, a probability is calculated for their colliding given the cross section and the fact that the sampling covers radians and one unit of . Given that the two radii are equal and that the other coordinates are irrelevant given the symmetries, the effective interaction range is zero. Furthermore, since any correlations from collisions or resonant decays are spread out over a wide range of and , this treatment should come extremely close to a true Boltzmann description even though particles are represented on a onetoone basis.
The algorithm is optimized by storing the information for each particle in a list ordered by radius. A second ordered list stores the list of pairs that will cross and is ordered by crossing times. The crossings are executed in order of time, and since the list is ordered, new crossings need only be calculated for the nearest neighbors of those particles that have crossed. When two particles cross, one needs to calculate the probability that they collide, or merge to form a resonance. This is related to the ratio of the cross section to the area of the cylinder over which the particles are spread, (assuming one unit of is being modeled). The exact probability is complicated by the relative angles of the particles and relativistic effects, and is given by:
(16)  
Here refers to the radial components.
All collisions before 25 fm/ are simulated, with decays being performed until they are exhausted. Weak decays are allowed to take place, except for charged kaons, pions and mesons. The point at which each particle had its last collision is recorded to be used for calculating spectra and femtoscopic correlations. A sampling of phase space points is displayed in Fig. 3 for particles with transverse momentum MeV/, and reveals a modest positive correlation between the outward position and time. This is opposite what one expects for an inwardly burning source whose emitting surface would move inward with time. Instead, it suggests an exploding source.
A glance at Fig. 3 reveals a strikingly different source than that extracted from Blastwave models. First, both the average lifetime and the duration of the emission are longer than those determined by blastwave models Retiere:2003kf (); Kisiel:2008ws (). Secondly, whereas the blastwave analysis of Kisiel:2008ws () suggested a negative correlation, there is a modest positive correlation in Fig. 3. Such a positive correlation was also seen in the AMPT model ampt (), which is based on a cascade picture for both partons and mesons. Since the hypersurface representing the transition from hydrodynamics to the Boltzmann approach has a negative correlation, this emphasizes the importance of accurately accounting for the breakup dynamics with a microscopic model. An underestimate of the emission duration for blastwave models is expected if the blastwave model employed an correlation of the wrong sign, as a positive correlation allows one to maintain a small ratio despite a longer emission duration. The mean emission time is associated with the ratio , where is the thermal velocity for longitudinal motion. Since blast wave analyses typically ignore shear effects and thus assume thermal motion is locally isotropic, they would overestimate if in fact the local momentum distribution is broader in the transverse plane than along the longitudinal direction. An overestimate of would lead to an underestimate of the lifetime.
The algorithm is both efficient and accurate. The procedure eliminates artifacts associated with particles interacting at a finite separation, and running on a single CPU, an event is performed in less than ten seconds. However, the approach has one numerical disadvantage. When calculating the value of at which two particles will reach the same radius, one must solve quartic equations. Numerical errors for such solutions are nonnegligible which occasionally lead to particles being propagated in such a way that violates the ordering by radius. This calculation is performed millions of times within a single event, and the violations tend to occur approximately once per every ten events. In such an instance, the event is abandoned. This does not seriously detract from the numerical efficiency, but such failures significantly complicated the construction of the code, as frequent error checking is required.
ii.3 Generating and Fitting Correlation Functions
Correlations can be generated via the Koonin equation,
(17)  
where is the pair’s momentum, is the spatial separation of the particles in the frame of the pair, and is the outgoing relative wave function. The probability of emitting a pion of momentum from spacetime point is , with being the coordinate in the pair frame. The source function is simply the normalized probability that two particles of the same momentum, , are separated by in the pair’s center of mass. The relative wave function incorporates quantum symmetrization, and both the Coulomb and strong interaction between the two pions. For the calculations shown here refers to one half the relative momentum as measured in the pair frame.
In practice, the Koonin equation is straightforward to implement. To calculate , one first extracts the subset of phasespace points from the output of the Boltzmann codes whose transverse momenta are within 5 MeV/ of . For every pair in the subset, one calculates for an array of values. The same set of pairs is used for every value of . The correlation function is then the average of for the pairs. Statistics for such calculations are greatly enhanced by the rotational and boost invariances, as every particle’s phase space points can be rotated and boosted so that it has zero longitudinal momentum and travels in a given azimuthal direction.
If one neglects the interpair Coulomb and strong interactions, the calculations can be greatly accelerated by calculating
(18) 
where the sum covers the particles in the subset used above, and is the position of the particle. For large , correlations can be generated by simply squaring the sum,
(19) 
Since one never has to evaluate a doublesum, this method is quicker than the alternate method. Although it can only be applied if one neglects strong and Coulomb forces, this method should be sufficient if one is generating correlations for the purpose of finding effective Gaussian source sizes. Both methods were tried for the calculation with the default parameters, with the comparison being illustrated in Fig. 4. Since the differences in the extracted Gaussian source sizes were small and the trends of interest are unlikely to be affected, the latter method was chosen. Although the calculations using the full wave functions are more realistic, it should be pointed out that experimental analyses have generated source radii by dividing the experimental correlation function by a dependent factor with the purpose of dividing away the effect of the Coulomb force in affecting correlation funcitons starhbt (); phenixhbt (). Since the Coulomb correction factors are based on isotropic Gaussian sources, the procedure is not exact. Now that the discrepancy between models and experiments are less than 10%, the errors introduced in this procedure should be reexamined. In particular, errors should be checked for radii at higher . For MeV/, source functions are highly anisotropic in the frame of the pair due to Lorentz contraction, with approaching five times , whereas the correction factors are built assuming that the source shape is isotropic in that frame.
To generate Gaussian radii, correlations were calculated on a threedimensional mesh in . These were then compared to predictions for for Gaussian sources. The radii were chosen to minimize the sum of the squared radii. Even though the correlations were remarkably nonGaussian for MeV/, the radii were remarkably robust, and did not change appreciably if one neglected the low points in the fit. Since the mesh and values were generated in the pair frame, the outward source size was then Lorentz contracted so that it represented the shape of the outgoing phase space density in the laboratory frame.
Iii Femtoscopic Ramifications of Adjusting the Equation of State, Viscosity and Initial Conditions
One of the prime motivations for interferometric measurements was the possibility of observing a longlived mixed phase, which could only happen if the equation of state was first order with a large latent heat. Some bag model parameterizations employed in the early days of the field had latent heats of several GeV/fm. If that were the case, extracted lifetimes at the AGS and SPS might have been several tens of fm/ depending on the amount of initial stopping Rischke:1996em (); Pratt:1986cc (). The longest lifetimes would occur for conditions where the interior energy density was initially at the peak value for the mixed phase. Since a mixed phase has zero sound velocity, there would no impetus for explosion, and instead the outside would emit hadrons like a burning log. Lattice calculations now preclude such equations of state, and indeed, no such long lived phases have been observed. Instead, in lattice calculations the speed of sound appears to dip down to about , before restiffening to at high temperatures Cheng:2007jq (). By including resonances in a hadron gas, the speed of sound is expected to be below . Given that initial energy densities at RHIC are well above those of the soft region, one should expect RHIC collisions to be more explosive and shorterlived that those at the AGS and SPS. Qualitatively, these expectations have been met. However, quantitatively describing the source sizes with full hydrodynamic models has proved elusive. Femtoscopy provides a sixdimensional test of any dynamical model, so it should not be surprising that reproducing experimental source sizes requires using a realistic equation of state, accurately modeling viscous effects and using correct initial conditions. We explore the impact of each of these three aspects of the modeling in the next three subsections. Results from the default parameter set are compared to results where an isolated parameter set has been adjusted. For each calculation, the initial energy density is adjusted so that the final phobosdndeta ().
Radial flow, and in turn spectra, are also affected by all the variations studied here. Table 1 presents the mean for pions, kaons and protons. Again, it should be emphasized that these calculations include all weak decays of hyperons and of the . To some extent, these decay products are subtracted from experimental analyses, which might lead to the model predictions under predicting the mean for pions. However, the calculations also neglect symmetrization effects on the pion spectra, which should lower the pions mean . Unfortunately, the uncertainties in the experimental values in Table 1 are rather large, mainly due to the fact that experiments measure in a finite range. The experimental values for mean do agree with the default calculation, within the large experimental uncertainties. A more meaningful comparison, which would involve comparing actual spectra in the measured regions, is outside the scope of this study, but Table 1 is certainly sufficient for evaluating the sensitivity of the spectra to the various model parameters studied here. It should be emphasized that the sensitivity of spectra to the equation of state has been considered by numerous authors, and that in Huovinen:2005gy (), sensitivity of the elliptic flow is also considered.
STARstarpt ()  

PHENIXphenixpt ()  
L=0  528  897  1310 
L=800 MeV/fm  433  714  1027 
L=1.6 GeV/fm  403  652  931 
406  659  945  
433  714  1027  
463  772  1116  
=0  408  664  957 
=2  433  714  1027 
=4  449  743  1081 
Initially isotropic  428  695  1012 
462  763  1107  
433  714  1027  
418  679  983  
CGC IC  447  741  1062 
Wounded Nucleon  433  714  1027 
Collision Scaling  482  806  1173 
iii.1 Adjusting the Equation of State
To study the sensitivity to the equation of state, we vary both the speed of sound and the width of the intermediate region, with the five equations of state being displayed in Fig. 5. For temperatures below 170 MeV, or equivalently for energy densities below MeV/fm, the equation of state is that of a resonance gas. For the intermediate region, , the equation of state has a constant speed of sound, . Above the intermediate region, the speed of sound is set to , to be consistent with lattice calculations at high temperature.
The equation of state can be softened by either increasing or decreasing the mixedphase value of . Figure 6 displays source sizes for the three values of : the default value of 800 MeV/fm, a soft value of 1.6 GeV/fm, and a hard value of zero. The default value was chosen to be crudely consistent with the behavior of lattice calculations which show a strong stiffening of the matter for energy densities rising from GeV/fm. The equation of state was also altered by adjusting the mixedphase value of from the default value of 0.1 to either a stiffer value of 0.2 or to a softer value of zero, which would correspond to a firstorder phase transition. The femtoscopic effect of varying the speed of sound is shown in Fig. 7.
As expected, softer equations of state lead to longer relative values of and relative to . Whereas the increase in signals an increase in the mean emission time, the increase in is indicative of a longer duration of the emission, or a more outsidein nature to the emission. The product of the three radii increases for softer equations of state. This is due to the increase in entropy associated with a softer equation of state (when compared at the same energy density). Although these variations in the equation of state are rather strong, doubling or , femtoscopic radii were affected on the level of 10%. Spectra are also affected by changes to the equation of state at the level of 10% as seen in Table 1. In particular softer equations of state lower collective radial flow which leads to lower values of the mean , especially for heavier particles.
iii.2 Adjusting Viscosities
Even modest viscosities signficantly modify the stress energy tensor. It has been proposed that shear viscosity can not fall below the KSS limit, kss (). According to Navier Stokes, at early times where the velocity gradient is , the KSS limit yields
(20)  
where the expressions involving assumed a free gas equation of state, . One expects for thermalization times near 1/2 fm/ where , that the correction to the longitudinal pressure is . If the viscosity is more than twice the KSS bound the value of can become negative. One expects a higher shear to accelerate the radial flow and result in lower values of and Paech:2006st (), as well as increased for heavier particles. These expectations have already been demonstrated by Romatschke romatschke ().
The default calculation presented here assumes that the shear viscosity is twice the KSS bound. This would yield negative at early times, if not for the mapping described in Eq. (6) which restricts the modification from shear to be less than the absolute value of the pressure. In the default calculation the initial condition for the stress energy tensor was set to this maximum value, with and , consistent with colorglass calculations Krasnitz:2002mn (). For pure noninteracting classical fields, the anisotropy would be even larger as is negative, and .
In addition to the default calculation, three modifications to the shear viscosity are illustrated here. First, the viscosity in the high density phase is set to zero. For the first modification the initial anisotropy to the stress energy tensor is also set to zero. For the second modification, the viscosity in the highdensity phase is doubled relative to the default calculation to four times the KSS bound. Due to the ceiling imposed on the viscous modifications, the higher viscosity only matters for times greater than 1 fm/. For the final modification, the default calculation is modified by setting the initial anisotropy of the stress energy tensor to zero. This mainly affects the expansion during the first one fm/. Since the IsraelStewart relaxation times tend to be fm/, memory of the initial condition is lost after that point.
The expectation for the femtoscopic radii are borne out by the results illustrated in Fig. 8. The default calculation differs from the zeroviscosity calculation by %. In particular, the ratio comes significantly closer to unity. It should be pointed out that the zeroviscosity calculation differs from many previous ideal hydrodynamic calculations in that the initial time was set to 0.1 fm/, whereas several other calculations used either 0.6 or 1.0 fm/, which would further increase the ratio. However, there is no physical justification for setting at early times, thus such an initial state seems unwarranted. This is discussed in more detail in Vredevoogd:2008id (). Modifying the initial anisotropy is similar in principal to altering the viscosity for early times.
Bulk viscosity is only expected to be significant near the critical region. In particular, the condensed fields may not be able to keep pace with rapidly changing equilibrium values. This can lead to a peak in the bulk viscosity in the intermediate energy region Paech:2006st (), which has been verified with analysis of lattice results Karsch:2007jc (). The divergence of the velocity , incorporates velocity gradients in all three directions is approximately one third (fm/) for fm/. For the default calculation, the peak value the in the intermediate region is . The magnitude of the effect is similar to what was derived in Paech:2006st (). For such a velocity gradient the trace of the stressenergy tensor is modified by a substantial fraction. Doubling the bulk viscosity can make the NavierStokes value of . For the calculations performed here, the mapping procedure of Eq. (6) saturates the size of the change in to be less than . However, when combined with shear effects individual components can fall below zero. It should be emphasized that the nonequilibrium effects that generate bulk viscosity, mainly nonequilibrium fields, may be very poorly represented by viscous formalisms. First, in the transition region, responses may be highly nonlinear, and secondly the field might not relax exponentially toward equilibrium as is assumed in IsraelStewart treatments. Thus, the study here can really only point at the qualitative impact of bulk viscosity on the dynamics, and ultimately on the femtoscopy.
The default calculation was modified twice, once by doubling the bulk viscosity and once be eliminating it. The difference of the three calculations, illustrated in Fig. 9, was modest. The effect on the stress energy tensor is similar to what one would get by softening the equation of state in the transition region. But unlike the changes to the equation of state investigated in the previous subsection, these changes do not affect the pressure at either high or low density. Hence, the impact on observables tended to be modest. Increasing the bulk viscosity only affected the femtoscopy at the level of a few percent. The bulk viscosity helped amplify the magnitude of the pulse in the energy density created near the soft region. However, the pulse largely dissipated later in the collision. The bulk viscosity had a slightly more significant impact on the mean as seen in Table 1 and in the initial energy density, which was adjusted to match , as seen in Table 2.
Bulk viscosity had a visible impact on the smoothness of the energy density profiles. Larger bulk viscosities appeared to lead to jagged and unstable profiles in the intermediate region. We speculate that this is driven by the fact that bulk viscosity effectively pushes the pressure vs. energy density to behave nonmonotonically, which could give regions where , which is the effective speed of sound squared, is zero or negative. It would be interesting to know whether such instabilities would appear in a more physically grounded treatment of nonequilibrium effects, such as one where the dynamics of nonequilibrium fields were treated in parallel to the hydrodynamic treatment by solving a coupled KleinGordon equation.
The significant sensitivity of the final femtoscopic source sizes to acceleration during the first one or two fm/ might seem surprising. The importance of early acceleration can be likened to an olympic sprint, where a head start of a few tenths of a second results in a difference of several meters at the end of the race. For this reason, it is imperative to understand the bulk properties, e.g., the stressenergy tensor, of matter even before thermalization.
iii.3 Adjusting Initial Conditions
For a rotationally and boost invariant calculation the choice of initial conditions involves choosing the initial transverse density, the initial stress energy tensor, and the initial time. For all calculations presented here, the initial time was chosen to be 0.1 fm/. Since the development of collective flow at such early times is driven by the initial stress energy, there is no reason to pick a later time, unless there were reason to expect and to be zero at early times. Since a little more than 0.1 fm/ is required for the nuclei to pass one another, and given that this is already an extremely short time relative to the overall expansion time, there is no motivation to pick an earlier time. Variations of the initial stress energy tensor, and in particular variations to the initial anisotropy of were considered in the previous subsection along with variations in the viscosity.
Three variations of the initial energy density profile were explored. The default calculation was that of the wounded nucleon model woundednucleon (). In this calculation the probability of a nucleon interacting is calculated as unity minus the probability it survives without interaction. The survival probability is calculated assuming the particle travels through the WoodsSaxon nuclear profile with a 40 mb cross section. For the thick part of the nucleus, this approaches participant number scaling. An alternative is collision scaling, where the energy density at a transverse coordinate is proportional to , where the thickness function is calculated by integrating the density of a nucleus over the coordinate. The third profile explored here is for the colorglass profile used in Drescher:2007cd (). In that case the energy density is chosen proportional to the minimum of and . For all three profiles, the thickness functions were found by convoluting two WoodsSaxon profile whose centers differed by an impact parameter of 2.21 fm, corresponding to the 5% most central collisions of Au+Au. The density profile was then averaged over the azimuthal angle to generate an approximate radial profile. For every calculation, the profiles were renormalized so that the resulting was , consistent with phobosdndeta ().
The collisionscaling profile was the most compact of the three attempted here, and resulted in the largest radial flow. This profile resulted in higher for protons and lower values of as seen in Fig. 10. The least compact profile was the default calculation which was based on the wounded nucleon model.
The initial energy density at the extreme periphery of the collision should be driven by collisional scaling since one can ignore the possibility of multiple collisions in that limit. None of the three profiles employed here obey this constraint, as the overall normalization is scaled in order to match the experimental value . If this constaint were replaced by the constraint that the density profile behaved correctly for small and , the experimental could instead additionally constrain the remaining parameters. The experimental has been used to argue that the profile appears more driven by participant scaling than collision scaling phobosdndeta (). However, these analyses are exceedingly simple and neglect many aspect of the expansion such as viscosity or longitudinal work. Since the temperature is not constant across the profile, the relation between entropy and energy densities is not linear, thus different profiles are reached if one believes the energy density should follow a given scaling vs. the entropy density, which is more related to the multiplicity. This underscores the importance of performing a global analysis of all observables, including elliptic flow, which is also affected by the choice of profile shape Drescher:2007cd ().
It is unfortunate that the initial energy density is not itself an observable. After adjusting the initial energy density for each parameter set to match the experimental , the final state observables tend to change at the 10% level or less for all variations studied here. However, the initial energy densities vary by more than a factor of two for these calculations as seen in Table 2. This variation has little to do with the variation of the asymptotic transverse energy for fixed multiplicity, which is also a % effect. The variation largely derives from differences in the longitudinal work in the expansion. The work is proportional to both and the time over which the system expands. For large viscosities, softer equations of state, or for initial conditions where is small, the longitudinal work is reduced, thus requiring smaller initial energy densities to attain a given final condition. For shear anisotropies and are also enhanced, which accelerates the expansion. This results in a reduced time for the reaction, which also reduces the longitudinal work. Changing the shape of the initial profile also changes the energy density, mainly for the trivial reason that a more compact initial profile requires a higher energy density to produce the same .
width of soft region in EoS  (GeV/fm) 
L=0  150 
L=800 MeV*  114.5 
L=1.6 GeV  104.5 
stiffness of soft region in EoS  
107  
*  114.5 
124.5  
shear viscosity in parton phase  
=0  289 
=2*  114.5 
=4  106.5 
initially isotropic init. cond.  148 
max. bulk viscosity in soft region  
=0  124 
=2*  114.5 
=4  109 
initial density profile  
CGC IC  136 
Wounded Nucleon*  114.5 
Collision Scaling  180 
Iv Summary and Outlook
The principal conclusion from these investigations is that the femtoscopic data from RHIC can be reproduced to within 10% with models combining viscous hydrodynamics and hadronic cascades. In particular, the ratio can be brought down close to unity. The failure of previous models appears to derive mainly from three shortcomings, all of which are related to underpredicting the explosivity of the collision. First, the equations of state were often too soft, using a firstorder phase transition. A stiffer equation of state is more explosive, and can lower the ratio. Secondly, previous treatments ignored acceleration before the thermalization time. From general arguments involving conservation of energy and momentum in the equations of motion of the stressenergy tensor, it should be clear that thermalization is not required for acceleration. In fact, longitudinal classical fields, which are far from equilibrium by definition, result in strong transverse acceleration. Finally, the previous treatments were based on ideal hydrodynamics. The effects of shear, as already demonstrated in romatschke (), increase the transverse pressure relative to the longitudinal pressure at early times, which of all the variations considered here, appears to be the most important. Bulk effects, were manifest in the final mean , but made remarkably little difference in femtoscopic radii. Previous treatments overpredicted the ratio by 40% or more Soff:2000eh (), a result confirmed if we run this model with a softer equation of state, without viscosity, and delaying transverse acceleration for the first fm/.
It would appear that improving models in all three areas mentioned above is required for rectifying the shortcomings. The default calculation, which includes all three such effects, provides a reasonable description to the data. Without including longitudinal acceleration, which requires a threedimensional model, it is unreasonable to expect better agreement and it is probably not meaningful to try to better fit the data by adjusting parameters. An additional area of uncertainty documented here comes from the choice of the initial profile, as a more compact source results in a more explosive source. One could reduce any of the three effects mentioned in the previous paragraph, then compensate for them by adjusting the initial density profile.
A second impression generated by this investigation is that it appears impossible to disentangle the various uncertainties mentioned above by focusing only on twopion interferometry. Spectra are sensitive to the same model features studied here, as evidenced by the mean values listed in Table 1. Elliptic flow observables, which require a higherdimensionality model than used here, can also be used to assist in understanding the collision. Hopefully, different observables will be relatively more sensitive to different facets of the model. Then by performing a coordinated analysis of numerous classes of observables, one should be better able to answer specific question and determine specific parameters. These analyses should also incorporate a greater set of femtoscopic measurements, which we list here:

Femtoscopy using particles other than pions. Heavier particles are more sensitive to collective flow due to their lower thermal velocity. Correlations between a heavy particle and a light particle, e.g. , are especially sensitive to the patterns of collective flow.

Noncentral collisions and collisions at different energies. Measurements have already been made as a function of the direction of the pair’s momentum relative to the reaction plane starhbt (). This information is rich in detail, but the meaning of the information is not yet understood. Additionally, there exists data at different beam energies. By studying the response to changes in the initial energy density, without changing the size, one should gain some leverage for disentangling some of the issues mentioned above.
Finally, it should be emphasized that, although the large discrepancy with the ratio has been eliminated, none of the variations studied here provided a completely satisfactory reproduction of the dependence of source dimensions. The data showed a modest fall of the ratio for higher , which combined with the constraint that the ratio must be unity for , gives a nonmonotonic behavior. Although the rise and fall are only of the order of 10%, the model calculations all showed monotonic behavior with . Furthermore, the model calculations tend to overpredict at low . This might be due to the assumption of boost invariance, which if relaxed, should provide corrections in the direction of the data. Finally, conspicuous by its absence, has been a comparison of the factors, which represent the fraction of pairs that are correlated. The model calculations overpredicted these factors, but without a better understanding of experimental details about acceptance of weak decay products and particle misidentification fractions, one can not as yet draw any conclusions.
Acknowledgments
Support was provided by the U.S. Department of Energy, Grant No. DEFG0203ER41259.
References
 (1) M. A. Lisa, S. Pratt, R. Soltz and U. Wiedemann, Ann. Rev. Nucl. Part. Sci. 55, 357 (2005).
 (2) A. Kisiel, W. Broniowski, M. Chojnacki and W. Florkowski, arXiv:0808.3363 [nuclth].
 (3) A. Kisiel, W. Florkowski and W. Broniowski, Phys. Rev. C 73, 064902 (2006).
 (4) F. Retiere and M.A. Lisa, Phys. Rev. C70, 044907 (2004).
 (5) J. Helgesson, T. Csorgo, M. Asakawa and B. Lorstad, Phys. Rev. C 56, 2626 (1997).
 (6) S. Pal and S. Pratt, Phys. Lett. B 578, 310 (2004).
 (7) P. Romatschke and U. Romatschke, arXiv:0706.1522 [nuclth] (2007); P. Romatschke, Eur. Phys. J. C 52, 203 (2007).
 (8) S. Pratt and K. Paech, Proceedings of the 22nd Winter Workshop on Nuclear Dynamics, La Jolla, CA, 2006. Ed. W. Bauer, R. Bellwied and S. Panitkin, EP Systema, Budapest (2006). [arXiv:nuclth/0604007].
 (9) M. Gyulassy, Yu. M. Sinyukov, I. Karpenko and A. V. Nazarenko, Braz. J. Phys. 37, 1031 (2007).
 (10) Q. Li, M. Bleicher and H. Stocker, Phys. Lett. B 659, 525 (2008).
 (11) U. W. Heinz and S. M. H. Wong, Phys. Rev. C 66, 014907 (2002).
 (12) W. Broniowski, M. Chojnacki, W. Florkowski and A. Kisiel, Phys. Rev. Lett. 101, 022301 (2008).
 (13) S. Pratt, Phys. Rev. C 77, 024910 (2008).
 (14) S. Weinberg, Gravitation and Cosmology, John Wiley and Sons, Inc. (1972).
 (15) W. Israel, Ann. Phys. 100, 310 (1976); W. Israel and J.M. Sewart, Ann. Phys. 118, 341 (1979).
 (16) A. Muronga, arXiv:0710.3280 [nuclth] (2007); A. Muronga, arXiv:0710.3277 [nuclth] (2007).
 (17) T. Koide, Phys. Rev. E 75 060103, (2007); Phys. Rev. E 72, 026135 (2005).
 (18) R. Baier, P. Romatschke and U. A. Wiedemann, Phys. Rev. C 73, 064903 (2006).
 (19) U. W. Heinz, H. Song and A. K. Chaudhuri, Phys. Rev. C 73, 034904 (2006).
 (20) H. Song and U. W. Heinz, Phys. Lett. B 658, 279 (2008).
 (21) D. Jou, J. CasasVázquez and G. Lebon, Rep. Prog. Phys. 51, 1105 (1988).
 (22) S. Pratt, Phys. Rev. C 75, 024907 (2007).
 (23) P.Kovtun, D.T. Son and A.O. Starinets, Phys. Rev. Lett. 94, 111601 (2005).
 (24) K. Paech and S. Pratt, Phys. Rev. C 74, 014901 (2006).
 (25) F. Karsch, D. Kharzeev and K. Tuchin, Phys. Lett. B 663, 217 (2008).
 (26) S. Cheng et al., Phys. Rev. C 65, 024901 (2002).
 (27) Z.W. Lin, C.M. Ko and S. Pal, Phys. Rev. Lett. 89, 152301 (2002).
 (28) J. Adams et al., Phys. Rev. C71, 044906 (2005).
 (29) K. Adcox et al., Phys. Rev. Lett. 88, 192302 (2002).
 (30) S. Soff, S. A. Bass and A. Dumitru, Phys. Rev. Lett. 86, 3981 (2001).
 (31) U. W. Heinz and P. F. Kolb, J. Phys. G 30, S1229 (2004).
 (32) D. H. Rischke and M. Gyulassy, Nucl. Phys. A 608, 479 (1996).
 (33) S. Pratt, Phys. Rev. D 33, 1314 (1986).
 (34) M. Cheng et al., Phys. Rev. D 77, 014511 (2008).
 (35) B.B. Back et al., Phys. Rev. C65, 061901R (2002).
 (36) P. Huovinen, Nucl. Phys. A 761, 296 (2005).
 (37) J. Adams et al., Phys. Rev. Lett. 92, 112301 (2004).
 (38) S.S. Adler, et al., Phys. Rev. C69, 034909 (2004).
 (39) A. Krasnitz, Y. Nara and R. Venugopalan, Nucl. Phys. A 717, 268 (2003).
 (40) J. Vredevoogd and S. Pratt, arXiv:0810.4325 [nuclth] (2008).
 (41) P.F. Kolb, J. Sollfrank and U. Heinz, Phys. Lett. B459, 667 (1999); Phys. Rev. C62, 054909 (2000).
 (42) H. J. Drescher, A. Dumitru, C. Gombeaud and J. Y. Ollitrault, Phys. Rev. C 76, 024905 (2007).