\thechapter Heavy Ion Collisions

# Relativistic Viscous Hydrodynamics for High Energy Heavy Ion Collisions

## Abstract

It has been over a decade since the first experimental data from gold nuclei collisions at the Relativistic Heavy Ion Collider suggested hydrodynamic behavior. While early ideal hydrodynamical models were surprisingly accurate in their predictions, they ignored that the large longitudinal velocity gradient meant that even small shear viscosity would produce large corrections to the transverse dynamics. In addition, much less was known about the equation of state predicted by lattice calculations of quantum chromodynamics, which predicts a soft region as the degrees of freedom change from quarks to hadrons but no first-order phase transition. Furthermore, the effects of late, dilute stage rescattering were handled within the hydrodynamic framework to temperatures where local kinetic equilibrium is difficult to justify. This dissertation presents a three-dimensional viscous hydrodynamics code with a realistic equation of state coupled consistently to a hadron resonance gas calculation. The code presented here is capable of making significant comparisons to experimental data as part of an effort to learn about the structure of experimental constraints on the microscopic interactions of dense, hot quark matter.

\field

Physics

\maketitlepage{dedication}

To my lovely wife and consistent source of inspiration, Sara.

\TOC\LOF
{doublespace}

## Chapter \thechapter Heavy Ion Collisions

It has been understood for over half a century that the protons and neutrons that make up atomic nuclei have substructure. The symmetry structure of the many new particles that were being discovered in the 1950’s and 1960’s led to an elegant explanation in terms of constituent particles that Gell-Mann coined as quarks. The quark model successfully predicted the Omega baryon, so named for being the only undiscovered combination of the three quarks (up, down, and strange) that made up all known particles at the time. Since that time, heavier quarks have been theorized and discovered. Notably, the top quark was discovered by the Tevatron collider at Fermilab in collisions of protons with anti-protons at extremely high relative momentum. Such high momentum collisions, of single hadrons or of nuclei, form the experimental basis of our knowledge about the interactions of quarks and the corresponding force carriers, gluons. The accepted theory for explaining these interactions is called Quantum ChromoDynamics (QCD), which has been studied in great detail and strenuously tested over the last half century.

Among the unique features of QCD is that each quark is confined by the strong force and can never be completely separated from at least one or two partner quarks. Confinement is borne out experimentally: despite the extreme amount of energy available to separate quarks in heavy ion collisions, collider experiments only observe mesons and baryons, made up of two and three constituent quarks respectively, and never a bare quark. In QCD, this is explained by the introduction of color charge. Isolating objects with non-zero color charge is forbidden, and since quarks carry a single color charge, they can never be observed in isolation. In the opposite regime, when quarks are extremely close, they interact only weakly, a result known as asymptotic freedom. This is exactly the opposite of Quantum Electrodynamics where the bare charge is infinite and renormalization is required to reproduce couplings observable at large distances as terms are added to the perturbation series.

Since the strong force increases dramatically as quarks become separated, there is the possibility of creating a new and interesting state of matter when the average distance between quarks becomes small. Creating this new state of matter, often called the Quark-Gluon Plasma (QGP), in which the quarks are freed from hadrons and form a plasma was one of the goals of the heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC), and the associated STAR (Solenoid Tracker At RHIC) and PHENIX (Pioneering High Energy Nuclear Interaction eXperiment) experiments [10]. While such a phase transition had been expected, clean experimental signatures that the transition has occurred have been somewhat difficult to isolate and, for instance, there is at present no evidence for behavior associated with a first-order phase transition. At this time, this is true even in the data from new RHIC initiative searching for the critical point in the QCD phase diagram at larger baryon chemical potential [11]. Still, significant experimental evidence suggests that the matter created at RHIC interacts strongly and collectively, notably away-side jet quenching [12, 1, 13, 14] and flow observables [15, 16, 17, 18].

Jets are formed in hadronic collisions at high relative momentum when there is a hard scattering of a pair of quarks or gluons. If the pair scatter away from the hadronic matter, the energy due to their separation grows rapidly. The separation energy becomes much larger than the energy required to generate additional quarks and so each quark or gluon begins to form quark-antiquark pairs until they are able to form a collection of colorless combinations. This process results in a group of high momentum hadrons at small relative momentum known as a jet. In proton collisions, jets are almost always found in back-to-back pairs with similar total energy due to momentum conservation. This is not the case in heavy ion collisions, as Figure 1 demonstrates for Au+Au collisions at a center of mass energy per nucleon pair () of 200 GeV. The figure shows the distribution of charged particles whenever there is a very high transverse momentum particle present shown as a function of the angular difference in the transverse plane compared to that particle. To assist in interpreting this plot, we define the coordinate system of heavy ion physics. The transverse direction is radially outward from the interaction point in the plane orthogonal to the motion of the colliding nuclei, where the direction along the the motion of the colliding nuclei is referred to as longitudinal. We will frequently refer to particles by their transverse momentum, , and their longitudinal rapidity, . The final coordinate is the azimuthal angle in the transverse plane, , meaning that the coordinate system in particle momentum space is cylindrical.

In Figure 1, the highest transverse momentum particle is defined to be at and only events with such a high transverse momentum particle are shown. Since events are selected on high transverse momentum, this particle is called the trigger particle and the associated jet is called the trigger jet. We expect a cluster of particles distributed around associated with the near-side jet, and another peak around corresponding to the away-side jet. This is strictly due to momentum conservation in the transverse plane. The particles in the trigger jet in the peak around zero are only slightly altered depending on the system size, but the away-side jet that appears as a peak at disappears for central Au+Au collisions, those that produce the largest total number of low momentum particles. This implies that the matter created in a heavy ion collisions is opaque to high momentum particles, though the theory of momentum transfer from high momentum quarks to a thermal QCD medium is not well understood [19].

The matter created in high energy heavy ion collisions also appears to interact long enough to cool. Figure 2 shows the chemical and kinetic temperatures indicated by blast wave fits to the transverse momentum dependence of the particle spectra as a function of the charged particle multiplicity. In its simplest form, the blast wave model assumes that all particles were emitted from a sphere of uniform density characterized by five parameters:

• A single kinetic temperature (),

• a single chemical temperature (),

• a linear velocity profile with unknown average velocity (),

• a spherical radius (R), and

• a baryon chemical potential ().

These five parameters are sufficient to predict the probability of observing a particular particle at a particular momentum, and therefore, they can all be determined by a fit to the spectrum low mass particles (pion, kaon, and proton). While this model is overly simplistic and does not include important physics, the parameters are a useful way of discussing trends in heavy ion collisions.

The chemical temperature extracted via the blast wave fit does not change as the number of charged particles or the beam energy increases and is consistent with the ratios observed in high energy proton-proton collisions. This trend does not carry over to the kinetic temperature indicated by the momentum distribution, which decreases as particle production increases. This softening of the spectra should not be associated with resonance production, since pions below 500 MeV/c were excluded, but instead should be related to the increasing importance of rescattering at increasingly temperatures. In peripheral collisions, which are shown in the figure as those with smaller multiplicities for the same beam energy, the spectra show no characteristic difference from proton-proton collisions in either the chemistry or the kinetics. This suggests that multiple scattering plays essentially the same role in these collisions. However, the effective kinetic temperature decreases smoothly for all beam energies as the hot region grows and the number of soft collisions, and therefore the system size, increase. This suggests that rescattering remains important for a considerable period after hadrons are created. A hydrodynamic model of heavy ion collisions must address the temperature range between kinetic and chemical freezeout. Modelers have taken two strategies – one can continue using a hydrodynamic model that takes into account the motion of the different particle species relative to one another; or one can couple the calculation to a gas model. In this work, we chose the latter.

The transverse velocity indicated by the blast-wave fit is further evidence of collective behavior in heavy ion collisions. Figure 3 shows the average velocity indicated by the blast-wave fit as a function of beam energy per nucleon pair for heavy ion collisions. A rapid increase in the average collective velocity is observed in the region around GeV and a continued steady increase toward the RHIC experiments. As a baseline, the STAR experiment finds that for proton collisions at GeV, the apparent average collective velocity is . The collective velocity in proton-proton collisions is not zero, as one might anticipate from a system expected to be too small to develop collective flow, but it is considerably smaller than as observed in heavy ion collisions [2]. This suggests that around a few GeV in beam energy the nature of heavy collisions changes, and collective response becomes critical to explaining the matter’s behavior.

In addition to the onset of significant outward flow, the data taken at RHIC show azimuthal anisotropy of flow as shown in Figure 4, which shows agreement between ideal hydrodynamic models and experimental data up to GeV. This anisotropic (or elliptic) flow is generated by the anisotropy in the initial state from the finite impact parameter. The elliptical interaction region means larger pressure gradients along the minor axis leading to more rapid hydrodynamic expansion. This additional flow produces additional particles at moderate momenta aligned with the minor axis across a large range of rapidities and therefore the second Fourier coefficient,

 v2=, (1)

increases where is the azimuthal angle relative to the impact parameter. The agreement between ideal hydrodynamic models was considered strong evidence for the formation of the Quark-Gluon Plasma [20, 21]. At the time, there were some caveats required. Ideal hydrodynamics overestimated the duration of the collisions as seen in the longitudinal source size [22], but systematic improvements to hydrodynamic calculations have resolved these discrepancies [23]. Elliptic flow is a very important observable in heavy ion collisions and will be a focus of Chapters Five and Six.

The truly novel feature of the Quark-Gluon Plasma is that quarks and gluons move freely and not in the bound particle states of a hadronic gas. While Figure 4 demonstrates that the matter created in a heavy ion collision exhibits hydrodynamic behavior, Figure 5 further shows that the anisotropic flow scales with the number of constituent quarks in the baryons and mesons observed in the final state. Scaling is not evident as a function of the transverse momentum but instead the transverse kinetic energy, , which corresponds to the kinetic energy available after creating the particle itself. If the anisotropic flow were generated in the liquid phase, translated to the hadronic matter by recombination, and subsequent rescatterings in the gas phase did not spoil the symmetry, the quantity would be expected to be the same for all particle species. The observation that this occurs in heavy ion collisions suggests that flow is developed prior to recombination.

Experimental data, shown in Figure 6, demonstrate that the observed chemical potential in midrapidity particles is very small compared to the temperature at the highest RHIC energies. This means that there is essentially no preference in the particle ratios between baryons and mesons other than that attributable to their differing masses. Therefore, we chose a hydrodynamic description of the evolution of the quark matter that ignores the effects of non-zero baryon number. In addition to reducing computational costs slightly, there is the theoretical issue that lattice calculations are not yet able to make reliable predictions for the equation of state of quark matter at finite chemical potential. Theoretical progress on this matter is ongoing [24, 25], but in light of the significant uncertainty we choose instead to limit our experimental observables to those at midrapidity to avoid these issues.

At the onset of this project, the state of the art in hydrodynamical modeling were viscous codes with trivial longitudinal expansion based on boost invariance [26, 27, 28]. Boost invariance is based on the observation that, in the center of momentum frame, small boosts should not influence observables since the original nucleons are at much larger rapidity than the created particles [29]. If this symmetry is observed exactly, then the system will have no variation in the coordinate

 η=12log(t+zt−z), (2)

where is time and is longitudinal distance relative to the symmetry axis, which we refer to as spatial rapidity. This is due to the correspondence of spatial rapidity to momentum rapidity,

 y=12log(E+pzE−pz), (3)

where is the longitudinal momentum. This correspondence is clearest when considering the product of a particle’s momentum with the collective velocity, as will arise in the phase space density,

 pμuμ=√m2+p⊥cosh(y−η), (4)

meaning that there will be a penalty for emission at a different rapidity that scales as . If one defines the proper time in a frame with constant as a replacement for the usual lab frame time coordinate, which is defined , one maintains an orthogonal set of coordinates with metric . Furthermore, the hydrodynamic equations of motion respect boost invariance, meaning that a system born into boost invariance maintains this symmetry at all times even in the presence of shear or bulk viscosity. In the limit of infinite beam energy, the symmetry is exact.

As noted earlier, for collisions at RHIC, boost invariance has proved a useful approximation for calculating midrapidity observables such as elliptic flow and particle spectra [26, 27, 28]. That said, the distribution of pions is still well described by a Gaussian in longitudinal rapidity. This was also the case at lower beam energies, though width increases to about 2.3 units of rapidity by GeV, as shown in Figure 7. While this means that boost invariance is increasingly valid, the range of validity was unknown as no three-dimensional simulation that included viscosity had been completed. In light of this, one goal of this project was to explore the validity of the boost invariant assumption in viscous hydrodynamic models regarding the prediction of midrapidity observables.

Taken in aggregate, the experimental evidence from relativistic heavy ion collisions points to the formation of a new liquid phase. This liquid remains in thermal contact long enough to develop significant radial and elliptic flow. However, the rapid longitudinal expansion means that even if the shear viscosity in the liquid phase is small, the viscous corrections to the transverse expansion will be important. For these reasons, we choose to model the high temperature regions with relativistic viscous hydrodynamics. Since larger systems cool further following the formation of hadrons, the lower temperature regions will be modeled as a hadronic gas. With this two-component model, we will investigate the sensitivity of our predictions at midrapidity to the addition of non-trivial longitudinal dynamics and uncertainties in the initial condition to hydrodynamics.

## Chapter \thechapter Theory of Fluid Dynamics

The experimental results discussed in Chapter 1 suggest that the matter created in heavy ion collisions can be described by hydrodynamics. However, thermal contact is maintained to lower temperatures than could be described by hydrodynamics. Because of this, our model will couple an hadronic gas at low temperatures to viscous hydrodynamics at higher temperatures. The focus of this chapter will be on hydrodynamic theory underlying the higher temperature model.

In this chapter, we pursue one of the possible methods of developing fluid dynamics. We consider the macroscopic, collective behavior of systems of free moving particles interacting with one another – a highly interactive limit of a gas. This will be of particular use when coupling the hydrodynamic phase to the gaseous phase where several results from this derivation will be of direct use. While this is a useful way of deriving the equations of motion of hydrodynamics, the resulting equations are more general than a derivation from kinetic theory might imply. Most notably, the quark-gluon plasma is expected to be strongly-interacting, while kinetic theory is weakly-interacting. The equations of motion apply to systems where entropy is conserved locally, collective motion is sufficient to describe the dynamics of interest, and the acceleration of that collective velocity comes from the absence or presence of material nearby. The equations are useful only when one can express an equation of state to close the set of equations. For the case of ideal hydrodynamics with no charges, this is to provide the pressure as a function of energy density.

### 1 Kinetic Theory and Fluids

In the case of relativistic heavy ion collisions, the underlying microscopic theory at high temperature should be Quantum ChromoDynamics for a bulk system. The present state of the theory has not produced a dynamical theory evolving microscopic degrees of freedom that would be needed to describe the observations from heavy ion collisions, but it has produced an equation of state for hot quark matter from lattice models (see Chapter 3) [6, 30].

The general structure of hydrodynamic theory can be derived from breaking the system into a system of fluid elements, each of which contains many colliding particles [31, 32, 33, 34]. If the particles collide relatively infrequently, the system is dominated by particle diffusion; whereas fluid dynamics considers a system for which particles collide on a shorter time scale than the system’s shape changes. The frequent collisions will tend to move the distribution of particle momenta toward the equilibrium distribution as collisions randomize each particles’ momentum. In the case that quantum corrections to this distribution are not important, this momentum distribution is the Maxwell-Boltzmann distribution, given by

 f(pμuμ)∝d6Nd3rd3p∝exp(−pμuμ/T), (5)

where is the phase space distribution, is the particle momentum, is the fluid velocity, and is the temperature.

The phase space distribution leads to the number density current and to the stress energy tensor. For the number density, one integrates over all momentum modes, which we will denote by

 ∫dω=∫d3p(2π)3p0, (6)

where is the relativistic energy. Using this notation, one can then define the number density current to be

 nμ=∫dω pμf(p⋅u), (7)

and the stress energy tensor to be

 Tμν=∫dω pμpνf(p⋅u). (8)

For a collisionless system, the phase space density is conserved in the frame that moves with the particle’s momentum, . If the particles are allowed to collide, these can be included as a source term to the equation called the collision integral. This is summarized by

 pμ∂μf(uμpμ)=C, (9)

where describes the effects of collisions. This equation is often referred to as the Boltzmann equation, though it has many names depending upon the addition of other effects. It is generally considered the starting point for kinetic theory and is fundamental to the description of gases.

As described above, hydrodynamics can be seen as a limit of kinetic theory where collisions do not allow the system to significantly deviate from the equilibrium distribution. If one then integrates the Boltzmann equation over momentum, one obtains

 ∫dω pμ∂μf=∂μnμ=0=∫dω C, (10)

where the expression is zero assuming that collisions are local and locally conserve particle number in the aggregate. If one takes the first moment of Boltzmann equation with respect to momentum, one obtains

 ∫dω pνpμ∂μf=∂μTμν=0=∫dω pνC, (11)

which is zero if the collisions also conserve momentum and energy. Higher moments of the Boltzmann equation also vanish exactly if the system never deviates from local equilibrium. This defines the stress energy tensor () from the perspective of kinetic theory which will serve as our connection to hydrodynamic theory, most notably when coupling the hydrodynamic module to the gaseous module.

In order to calculate the connection with hydrodynamics, we consider a small volume of particles distributed according to the equilibrium distribution with no collective motion. If we calculate from kinetic theory as in Eq. 8, we obtain

 Ttt(→r)=∫d3pp0p20f(→p,→r)=ϵ, (12)

where is the energy density. Since this frame is defined by having zero net momentum, all off-diagonal elements of the tensor are zero at equilibrium in this frame due to symmetry. Spatial elements along the diagonal, however, are proportional to the average of a momentum squared, which is non-zero and given by

 Txx(→r)=∫d3pp0p2xf(→p,→r)=P, (13)

where is the pressure. The identification of this quantity as the pressure can be understood by noting that the linear momentum transferred to the boundary of the fluid element and the frequency of collisions with the boundary are both proportional to the momentum [31]. This will be confirmed later when we investigate the conservation of stress-energy in the following section.

One can also note directly from Eq. 12 and Eq. 13 that the equation of state for a massless gas is entirely determined. From relativistic kinematics one knows that

 p20−→p⋅→p=m2, (14)

where is the mass of the particle. But then if one examines the trace of ,

 Tμμ=Ttt−Txx−Tyy−Tzz=∫d3pp0(p20−p2x−p2y−p2z)f(→p,→r)=0, (15)

where the final step uses the fact that the particles are massless. Changing from the kinetic theory description to the language of hydrodynamics, Equation 15

 ϵ=3P, (16)

which suffices to provide an equation of state for ideal hydrodynamics. Note also that if the mass is non-zero, then the trace of the stress-energy tensor is positive and the pressure at constant energy density is reduced. This is relevant for temperatures larger than the pion mass where more mesonic and baryonic states become relevant and the equation of state begins to soften (Chapter 3).

### 2 Ideal Hydrodynamics

The equations of relativistic hydrodynamics can be determined by adding the requisite relativistic structure to the stress-energy tensor. The only available non-scalar quantities are the velocity of the frame in which the fluid element has no collective velocity and the metric tensor. No derivatives are permitted in the stress energy tensor since their inclusion would require the introduction an additional microscopic length scale over which we assume that the stress energy tensor will not vary. Easing of this assumption will come in the form of viscous corrections to be discussed in later sections. From our earlier symmetry consideration, the stress energy tensor has no off-diagonal structure and trace reads in the frame with no collective velocity. This is sufficient to determine the relativistic structure [32], which for ideal hydrodynamics is given by

 Tμν=uμuν(ϵ+P)−gμνP, (17)

where is the collective velocity and is the metric tensor. We take the metric to be in Minkowski space, and transforms as a vector and has unit measure, . Implicitly, we have chosen the Landau frame where the collective velocity moves with the energy density, and not the Eckhart frame where the collective velocity moves with particle number. This is due in part to interest in systems with zero net charge where particle number, and therefore the Eckhart frame, would not be well defined.

The equations of motion of ideal hydrodynamics are then just the conservation of the stress energy tensor as from Eq. 11. The portion of temporal component along the collective velocity is given by:

 0 = uν∂μTμν=uν∂μ[uμuν(ϵ+P)−gμνP], 0 = Dϵ+(ϵ+P)∂μuμ, (18)

where is the comoving time derivative. To obtain the familiar non-relativistic limit, assume that is small and that . This yields the familiar expression:

 ∂tϵ+→v⋅→∇ϵ=−(ϵ+P)→∇⋅→v. (19)

The divergence of the collective velocity is often called the expansion rate and can be related to the change in volume of the fluid element by . Equation 19 can then recast as which is the fundamental thermodynamic relation for a system with constant entropy. For this reason, the first equation of ideal hydrodynamics is often framed in terms of the local conservation of entropy. Furthermore, it underscores the more general applicability of hydrodynamics to statistical systems at or near the maximal entropy state.

The Euler equation is related to the conservation of the spatial elements of the stress-energy tensor. To extract this independently of the first equation we use the projector which selects out the spatial indices in the frame of the matter. Applying this to the equation of stress-energy conservation yields

 0 = Δαν∂μTμν=∂μTμα−uα(Dϵ+(ϵ+P)∂μuμ), 0 = (ϵ+P)Duα+uαDP−gαμ∂μP. (20)

When taking the same limit as before, one obtains

 (ϵ+P)∂t→v+(ϵ+P)[→v⋅→∇]→v=−→∇P−→v∂tP, (21)

which reduces to in the frame of the matter. This should be interpreted to mean that acceleration of the frame of the collective velocity is due to the pressure gradient observed in that frame.

### 3 Viscous Hydrodynamics

#### 3.1 Navier-Stokes Hydrodynamics

The equations of ideal hydrodynamics are derived from the assumption that local equilibrium is exact, or equivalently, assuming that the mean free path of the microscopic particles is identically zero. If this is the case, a uniform density fluid with a plane of static fluid adjacent to a plane of fluid in motion would be a stable condition; the two planes would never affect one another because only the divergence of the collective velocity, rather than a mixed partial derivative like , enters into the equations of motion. That is, pressure is a force orthogonal to the surface of the fluid element in the rest frame; shear viscosity allows force along the surface.

Experience with liquids indicates that this is not a physical conclusion; if one pulls a sheet of material through a fluid, nearby fluid begins to move with the sheet. Friction in the fluid allows linear momentum to be transferred between adjacent fluid elements and their collective velocities would tend to equalize. Navier-Stokes hydrodynamics introduces a viscous correction to the spatial elements of the stress-energy tensor from such a frictional force [31]. Physically, this should be separated into a traceless contribution and an effect on the trace, which correspond to shear () and bulk () viscosity respectively. The viscous components are proportional to the velocity gradient with a corresponding transport coefficient that characterizes the amount of internal friction in the fluid. In the traceless case, this the shear viscosity. These results can be summarized as

 δTij∝η∂vi∂xj,  δTij∝ζδij→∇⋅→v. (22)

These terms are collected into a correction to the stress energy tensor – – and referred to as the shear tensor. The shear tensor is mostly clear written in the fluid frame where it takes the form

 πij=−η(∂ivj+∂jvi−23δij→∇⋅→v), (23)

where if and is otherwise zero, and the latin indices are meant to indicate that this applies only to spatial coordinates. In the frame of the matter, the ideal part of the stress energy tensor is only the outward pressure, , and Eq. 23 is added as a correction to the equations of motion (Eqs. 19 and 21). These additions result in the Navier-Stokes equations of viscous hydrodynamics. Note that adding these corrections relativistically is non-trivial as one needs to correct for both the frame motion in both the derivatives and the velocities. This will be addressed fully when discussing the Israel-Stewart equations of motion.

#### 3.2 The Maxwell-Cattaneo Equation

Solving the Navier-Stokes equations is notoriously difficult as the differential equations are parabolic, for which solutions tend to be unstable, and famously unsolved in the general case. In addition to being difficult to solve, the relativistic Navier-Stokes equations are not necessarily causal as the system responds instantaneously to changes in the velocity gradients.

While it is not clear that either of these difficulties are important for heavy ion physics, both can be addressed within the Maxwell-Cattaneo framework [32]. Essentially this framework assumes that instead of defining the change to the effective pressure to be exactly proportional to the velocity gradients, it instead relaxes toward the Navier-Stokes value. This means that the shear tensor can no longer be computed directly from the other dynamical variables or their derivatives. Therefore, the shear tensor itself is promoted to a dynamical variable and must be tabulated separately. In addition, the shear tensor now must have its own equations of motion. The form of these equations of motion in the frame of the matter are taken to be of the form

 ∂t(πijση)=−(πij−πij(NS))σητπ, (24)

where is the Navier-Stokes shear tensor given by Eq. 23, is a new transport coefficient that characterizes the time scale on which the system approaches Navier-Stokes viscosity, and is a scaling factor required to guarantee the Second Law of Thermodynamics [35] (see the following subsection).

Historically, this equation informed the search for an extension of hydrodynamic theory that would would be hyperbolic and causal in the general case. The equations of motion for the shear tensor can be derived in this form from a gradient expansion in kinetic theory or from general entropy production considerations as will be shown in the next two subsections. That said, the new transport coefficients, and will not be independent of one another complicating the correspondence of the Israel-Stewart theory to the Maxwell-Cattaneo equation and its interpretation as a relaxation equation.

#### 3.3 Israel-Stewart Hydrodynamics and Entropy

Since the system is away from equilibrium, the entropy should be reduced from its equilibrium value. One expects the entropy should increase quadratically as the shear tensor increases simply from symmetry at equilibrium so the lowest order correction to the entropy current would be

 sα=uα(seq−βπμνπμν), (25)

where is a thermodynamic quantity to be determined [35]. The second term acts as an entropy source and so it must have positive four-divergence relative to the dynamical entropy from the Navier-Stokes continuity equation [32, 36]:

 0 ≤ 1Tπμν∇<μuν>−∂α[uαβπμνπμν], 0 ≤ πμν[1T∇<μuν>−2βDπμν−πμνDβ−πμνβ∂αuα], (26)

where is the traceless part on the velocity gradient and is defined as

 ∇<μuν>=∇μuν+∇νuμ−23Δμν∂αuα, (27)

with

 Δμν=gμν−uμuν,    ∇μ=Δμν∂ν. (28)

This traceless part of the velocity gradient is just the Navier-Stokes modification to the pressure except the factor of the shear viscosity: . The constraint in Eq. 26 must be guaranteed for any expansion with any physical transport coefficients, which can be achieved if one requires that

 πμν=ηT[1T∇<μuν>−2βDπμν−πμνDβ−πμνβ∂αuα]. (29)

This is an equation of motion for the shear tensor, where one should recall that is the time derivative in the frame of the matter. If is taken to be zero, the Navier-Stokes form for the shear tensor is recovered exactly. In addition, it can be transformed into a relaxation equation with the same structure as Eq. 24 via the appropriate choice for the new transport coefficient . To show this, one can rewrite Eq. 29 as follows

 Dπμν+πμνDβ2β+πμν2∂αuα=−(πμν−πμν(NS))2ηTβ, (30)

which leads to two equations for the two unknowns:

 β=2ηTτπ,     σηD(πμνση)=Dπμν+πμνDβ2β+πμν2∂αuα. (31)

Since the relaxation time is a transport coefficient to be calculated from the underlying microscopic theory, the left part of Equation 31 gives a value for . The right part of Equation 31 can be recast as a single, local conservation equation

 D(ση√sβ)ση√sβ=0, (32)

which is guaranteed if . Returning to Equation 25 for the viscous entropy current and noting that the probability for a given fluctuation of stress-energy is proportional to , then the variance of in a given volume V is

 V<π2xy>=12β=sσ2η2, (33)

meaning that behaves like the variance of .

#### 3.4 Israel-Stewart Hydrodynamics and Kinetic Theory

The derivation connecting Israel-Stewart hydrodynamics with kinetic theory is rather involved and the details are beyond the scope of this document, but a quick sketch is helpful for conceptual reasons. In addition, the form of the phase space distortion shown in this section will be useful later when connecting the hydrodynamic fireball to the surrounding hadronic gas.

In establishing the connection between kinetic theory and ideal hydrodynamics, we discussed the equilibrium phase space density and the Boltzmann equation. In doing so, we discussed the effect of collisions, which are a source (or sink) to the conservation of phase space density. We argued that if the collisions happened often enough, then the system would remain in equilibrium. This assumption along with energy-momentum conservation was enough to produce ideal hydrodynamics.

Further, when discussing the Navier-Stokes viscous corrections hydrodynamics, we mentioned that this can be thought of either as friction in the fluid from a macroscopic perspective, or as the effects of expansion between collisions. In particular, if one imagined a box of particles with an arbitrary distribution of initial momenta but zero total momentum, and then investigated the system at later times, the collisions would tend to relax toward the momentum distribution and the system would exponentially approach the equilibrium distribution [37, 33, 38, 39, 34]. This can be summarized as

 C=−pμuμf−f0τπ, (34)

where is the relaxation time scale.

Now, hydrodynamics will not generally apply to situations where the system is far from thermal equilibrium so we concern ourselves only to the case of small deviations of the phase space density. Using the fact that this deviation should disappear in equilibrium and should only be a function of the available degrees of freedom, the simplest ansatz was proposed by Grad [40] who found that the distortion to the phase space density should be

 f=f0[1+παβpαpβ2(ϵ+P)T2]. (35)

This form of the deviation and the equation of motion for the shear tensor come from a detailed analysis of higher moments of the Boltzmann equation [41, 33, 38].

The important revelation from this approach is the possibility of additional terms that do not produce any entropy [32]. Such terms would not appear in a derivation based solely on the form of the entropy current. These additional terms tend to require even more transport coefficients many of which do not have obvious physical descriptions. The exception is the coupling of the shear tensor to the vorticity tensor, which in the frame of the matter is

 Ωij=∂ivj−∂jvi, (36)

would arise naturally from hydrodynamic considerations and does not require additional transport coefficients. This term has been included in our calculations but our studies confirm others’ results that it is not an important driver of dynamics for smooth initial condition investigations of hydrodynamical flow.

At this point, we have discussed all the important physics related to viscous hydrodynamics for the purposes of producing a simulation of zero charge heavy ion collisions. The details of the equations of motion to be integrated have been reserved for Chapter 4.

## Chapter \thechapter Equation of State and Transport Coefficients

While hydrodynamics is a theory that applies to a wide variety of systems, both strongly and weakly coupled, the microscopic details of a fluid need to be included through the equation of state and the transport coefficients. For heavy ion collisions, two physical regimes are relevant as discussed at length in Chapter 1. At low temperature, the equation of state should be calculated as a non-interacting gas of mesons and baryons; while at high temperature, lattice QCD should be used. These predictions may not be compatible in the region near the transition between hydrodynamics and the hadron resonance gas, and care must be taken to ensure that the model is self-consistent at the boundary. This is especially important for the equation of state where, for instance, discontinuity in the pressure as a function of energy density would lead to hydrodynamic instability. While this temperature region is still under active investigation, recent lattice results suggest that the discrepancy between the two theories is not large and a composite description is possible. This is discussed in the first section of this chapter.

In the second section, we discuss the determination of transport coefficients in linear response theory. At this time, lattice calculations are not developed to the point where they can make predictions about transport coefficients using classic results like the Kubo relations. Because of this, our results will be from the weakly coupled theory at low temperature, the general expectation that shear viscosity should be small near the phase transition, and scaling arguments from scalar field theories of the high temperature phase.

### 4 Equation of State

#### 4.1 Hadron Resonance Gas

The equation of state of the hadron resonance gas model is developed as a gas. If one assumes that the gas is non-interacting, one can simply generate a partition function using all the available particle states, which are measured experimentally. This ignores the finite volume or mean field effects that are well known to be important for lower energy heavy ion physics but are expected to be less important corrections at energies where zero baryon number hydrodynamics apply. Data on particle states up to several GeV/ in mass are available from the Particle Data Group. We use resonances up to masses of GeV/ as the addition of higher mass resonances no longer affects the equation of state up to temperatures of MeV. This choice leads to a nice agreement with lattice results in the transition region [42].

The partition function is calculated from the usual product of available mass states weighted by a common temperature. The contribution to the log of the partition function for a mass state is given by

 1VlnZi=∓di2π2∫∞0dkk2ln(1∓zie−√k2+m2i/T) (37)

where is the spin degeneracy of the state and are the fugacities summed over all charges [6]. One can then use standard thermodynamic relations to obtain the equation of state. For instance the pressure and energy density are

 P=TV∑ilnZi,  ϵ=T2V∂∂T∑ilnZi. (38)

From these, the entropy density can be calculated using the thermodynamic identity, .

To briefly summarize the results of this procedure, the hadron resonance gas acts like a pion gas for temperatures less than 100 MeV. However, the increasing number of baryon and meson states play an increasingly important role as the temperature increases. This produces an order of magnitude increase in the trace anomaly (or interaction, ) by 150 MeV. In fact, inclusion of the resonances in the mass range provides a factor of 3 in the trace anamoly by MeV [6, 30]. From a hydrodynamic perspective, this yields a soft region – with lower speed of sound – as a heating system would expend part of its energy in changing the degrees of freedom instead of increasing the pressure. This behavior is a milder version of a latent heat which would be associated with a first order phase transition, where there would be no increase in pressure over a region of increasing energy density.

#### 4.2 Lattice QCD

In Lattice QCD, one also sets out to calculate a partition function as a function of temperature but within a field theory rather than a non-interacting gas. The details of this procedure are complex and involve computing path integrals over the gauge field on a discretized hypercube and extrapolating to the physical pion mass. The computational effort required for each time step is immense and computations are performed on large clusters designed specifically for the task. Developments in this area now allow for the computation of twelve time steps for spatial grid points [6, 30]. Different temperatures are investigated by altering the lattice spacing, and several spacings must be investigated to quantify convergence and protect against discretization effects.

While the lattice has produced some incredibly useful results for finite temperature QCD matter, it has important limitations. The “sign problem” complicates direct investigation of the equation of state away from zero chemical potential where there is the possibility of a critical point. There the phase transition might shift from a smooth cross-over seen at zero chemical potential to a first-order phase transition. Extrapolation toward this critical point from lattice data has been attempted via Taylor expansion [25, 24]. While one does not expect this expansion to be valid far from zero chemical potential, especially when approach a critical point, recent results suggest tentative agreement with other theories.

Furthermore, it is not possible to use a physical quark mass for the up and down (the lightest) quarks. The up and down quarks have physical masses of around 3-4 MeV and the strange quark around 200 MeV, but the typical lattice calculation take the mass of the up and down quarks to be a tenth that of the strange quark. This leads to a pion mass that is significantly too large. This significantly alters conclusions about the pressure in this region and generally makes interpretation of the lattice results significantly more difficult as one must extrapolate to the physical meson and baryon masses in order to extract corrected thermodynamic quantities. Despite of all this, as one finds in Fig. 8, there is good agreement with the Hadron Resonance Gas model for the temperature region in which one hopes both might be valid.

Also included in Fig. 8 is a parameterized fit to the unit-less interaction over the full temperature space. In this case [6], the fit function reproduces the HRG interaction within 7% at temperatures below 100 MeV, and calculation of the pressure via the integral method

 P(T)T4=∫T 0dT′T′I(T′)T′4 (39)

finds deviations of only 2% and other collaborations find comparable results [30].

#### 4.3 Cross-over Region

The general structure of our model is a hydrodynamic region defined by temperatures above some threshold surround by a gas of interacting hadrons. Practically, dynamic coupling of these two calculations is difficult as pressure fluctuations in the gas might easily result in extreme hydrodynamic instability. Instead, one runs the hydrodynamic simulation for a larger region of configuration space under the assumption that regions at much lower temperatures do not have undue influence on the hydrodynamical portion and that near the switching temperature the models are equivalent.

As a direct consequence, modelers should be fastidious in ensuring that the equation of state for the gas is used as precisely as possible for those regions. While lattice calculations provide a fit that includes data points from the HRG equation of state and the lattice, this may not be sufficient to ensure self-consistency. In fact, one finds that for MeV, lattice fits can produce energy densities that differ from the gas’ energy density by 10-20% without fundamentally different behavior in the interaction measure, . If left uncorrected, the hydrodynamic code would improperly emulate the gas at the boundary, and the model of particle production would fail to be self-consistent.

Since one might want to use an arbitrary lattice equation of state that might not reproduce the gas well enough, we perform a continuous merge between the two equations of state over a temperature range starting at the hadronization temperature () up to the temperature where only lattice results are used (). In particular, we wish to maintain a continuous speed of sound for hydrodynamic stability, which can be written in terms of the unit-less entropy density () as

 c2s=dPdϵ=13+TσdσdT (40)

which can be derived using standard Maxwell relations. This is a clarifying expression as it points out that entropy increasing faster than is associated with softness in the equation of state, but also means that is a useful variable for merging equations of state, as continuity and smoothness ensure the continuity of the speed of sound.

We perform the merge in the unit-less entropy density using a linear weight function that goes smoothly and continuously between the hadron resonance gas entropy () and the lattice entropy ():

 σ(T) = [1−w(T)]σH(T)+w(T)σL(T) (41) w(T) = 12[tanh(tan[π2(2TL−TTL−TH−1)])+1], TH

and is zero for and unity for and importantly also has zero slope at both and to maintain continuous speed of sound. Figure 9 shows the result of the procedure compared to the original equations of state and demonstrates that the procedure is continuous and smooth.

Like the hadronization temperature, the lattice-only temperature is a parameter of the model and should not be so large that gas data is being used at temperatures well above its applicability or so small that the entropy increase becomes too rapid. Figure 10 shows the speed of sound squared in the region of the merging procedure. Since the entropy of the hadron gas lags behind the lattice, the merging procedure gives a steeper increase for temperatures just greater than the hadronization temperature. This manifests as an even smaller speed of sound in this region. Lowering the lattice temperature while leaving the hadronization temperature the same would results in a more extreme dip in the speed of sound squared. Choosing a lattice temperature too close to the hadronization temperature could lead to an extremely soft region in the equation of state and cause hydrodynamic instabilities.

For completeness, we note that the entropy density is the temperature derivative of the pressure, , so knowing the pressure of the hadron gas at and the entropy density at all temperatures allows one to get the pressure at all temperatures. Then the energy density is found from . We calculate each quantity at an interval of 1 MeV in temperature and a bookmarked search is performed each time a thermodynamic quantity is needed. Since our model does not solve the Riemann fan and cannot handle shocks, we take great care to avoid interpolation errors in the equation of state. A cubic spline [43] is performed on each of the thermodynamic variables and the interpolated values from this are used including the speed of sound squared which must be calculated as part of the splining procedure. This procedure involves a minor amount of preprocessing but the bookmarked search algorithm is unchanged and a minor part of computational time.

### 5 Transport Coefficients

As discussed in the previous chapter, viscous hydrodynamics introduces the effects of non-equilibrium to the collective response of material. In the case of shear viscosity, this can be thought of as friction within the fluid or a finite mean free path, if the underlying theory is a classical gas. The form of the correction is a small change to spatial elements of the stress-energy tensor due to velocity gradients, for example,

 δTxy=η(∂xvy+∂yvx). (43)

This can be viewed in terms of Linear Response Theory where is the operator being varied, is the coefficient of response, and the velocity gradients are the small field.

For the variation of some operator , the variation of its expectation value due to the field can be expressed in terms of the variation of the wavefunction as

 Missing or unrecognized delimiter for \right (44)

where is a susceptibility and F a small external field. The change the wave function can be viewed as the result of some potential, which we assume to be related to the field by , leading to

 |δψ>=−iℏ0∫−∞dtV(t)|ψ0>=−iℏ0∫−∞dtB(t)F|ψ0> (45)

so that the variation of the original operator can be written

 δ⟨A⟩ = Missing or unrecognized delimiter for \left (46) χ = Missing or unrecognized delimiter for \left (47)

In the case of shear viscosity, we are interested in the response of stress-energy to the velocity gradient. The potential associated with this is

 V=∫d3rT0jri∂vj∂ri, (48)

where latin indices indicate spatial coordinates. This form for the potential can be viewed as the change to the energy due to boosting to the frame of the velocity gradient [44]. This means that Equation 47 gives the shear viscosity to be

 Missing or unrecognized delimiter for \left (49)

Inserting unity in the form and the integrating by parts with the correlation at minus infinity assumed to be zero,

 η = −iℏ0∫−∞tdt∫d3r⟨[Txy(0),∂tT0j(r,t)]⟩ri, (50) η = iℏ0∫−∞tdt∫d3r⟨[Txy(0),Txy(r,t)]⟩, (51)

where we have used the conservation of stress-energy, integrated by parts in the spatial coordinate disregarding the correlation at infinite range, and then noted that symmetry requires that some of the stress-energy tensor elements be the same. The commutator in Equation 51 can be eliminated:

 0∫−∞t dt⟨[Txy(0),Txy(t)]⟩ = 0∫−∞t dt(⟨Txy(0)Txy(t)⟩−⟨Txy(0)Txy(−t)⟩), (52) = ∞∫−∞t dt⟨Txy(0)Txy(t)⟩,

where we’ve ignored the spatial integration for clarity.

Producing the classical Kubo relation from this expression requires evaluating the thermal average, which is

 G(t) = ⟨Txy(0)Txy(t)⟩=Tr[e−βHTxy(0)Txy(t)], (53) = Tr[e−βH Txy(0) eiHt/ℏ Txy(0) e−iHt/ℏ].

After analytically continuation, for complex is symmetric about :

 G(z+iβℏ/2) = Tr[e−βH Txy(0) eiHz/ℏ−βH/2 Txy(0) e−iHz/ℏ+βH/2], (54) = Tr[e−βH Txy(0) e−iHz/ℏ−βH/2 Txy(0) eiHz/ℏ+βH/2]=G(−z+iβℏ/2),

where the cyclic symmetry of the trace is used to rotate around in Equation 54. Note that this result also means that . Together these results mean that integration along the positive real axis can be performed by a complex contour integral. The integration region is taken to be rectangular with corners at and extending to infinity in the positive real direction. The region is closed by a curve at infinity which disappears as long as correlations decay faster than for large . The contour integral in full is zero since the contour contains no poles, so

 0=∮dz (z−iβℏ/2)G(z), (55)

where the integrand is chosen to be odd about . This leaves only the integral along the positive real axis and the return integral along :

 0 = ∞∫0dt (t−iβℏ/2)G(t)+0∫∞dt (t+iβℏ/2)G(t+iβℏ), 0 = ∞∫0dt (t−iβℏ/2)G(t)−∞∫0dt (t+iβℏ/2)G(−t), 0 = ∞∫0dt t[G(t)−G(−t)]−iβℏ2∞∫0dt [G(t)+G(−t)]. (56)

This means that the shear viscosity can be written in terms of an anticommutator instead of a commutator,

 η=β2∞∫0d4r⟨{Txy(0),Txy(r)}⟩. (57)

In the classical limit, the stress-energy operators commute and Equation 57 becomes

 η=β∞∫0d4r⟨Txy(0)Txy(r)⟩, (58)

which is the familiar Kubo relation. Larger values of are then due to the persistence of correlations due to fluctuations and vice versa. This is especially clear in the case of a massless gas that loses correlation exponentially due to a finite time between collisions. Using Equation 13 for the stress-energy tensor of a massless gas, the correlation of stress-energy is

 ∫d3x⟨Txy(0)Txy(x)⟩=∫d3p(2π)3(pxpyp0)2e−t/tc−βp0, (59)

where is the average time between collisions. This integral and the time integral in the Kubo relation can be performed analytically and the result is a relationship between the collision time and the shear viscosity for a massless gas given by where P is the pressure [45]. The direct relationship between shear viscosity and relaxation time confuses the interpretation of the Israel-Stewart equation for viscous corrections as a relaxation equation. Regions where the shear viscosity is large are also regions where stress-energy is slow to respond to changes in velocity gradients. Other theories predict different relationships between the shear viscosity and the relaxation time, for instance in strongly coupled symmetric Yang-Mills theory the relaxation time is a factor of two smaller with the same temperature dependence.

In principle, Equation 58 could be evaluated on the lattice to give the shear viscosity of the quark matter for all relevant temperatures in the hydrodynamic calculation. Currently, this is impossible due to numerical noise though such a calculation remains an possibility given even more computational resources. For lower temperatures, it is possible to use kinetic theory to estimate the shear viscosity. Viewed in terms of the unitless ratio to the entropy density, , this tends to rise very quickly with decreasing temperatures as is typical for gases below their critical temperature – for instance, jumping from unity at around T = 100 MeV to five by T = 60 MeV for a massless pion gas [46]. Such temperatures are by design outside of the hot region and tend to have a small effect on the output from the hydrodynamic module. We investigated the effect of changing the behavior of shear viscosity at temperatures below the freezeout temperature combined with a constant shear viscosity to entropy density ratio at high temperature, . We considered three low temperature scenarios: constant , constant below T=130 MeV, and increasing as a Fermi function up to 0.6 at low temperature. This is still less than the value predicted by the hadron resonance gas model but still runs stably for reasonable parameter sets. None of these runs showed significant differences in the produced elliptic flow which justifies the choice to exclude this correction from the set of first parameters to investigate. Other studies in this area have found that low temperature shear viscosity can have an effect if it is allowed to rise for temperatures slightly above the freezeout temperature [47] most likely due to phase space effects during particle generation as others have found that high temperature shear viscosity scaling matters much more [48].

In the high temperature region, we also expect the shear viscosity to entropy density ratio to rise. A scaling argument for this temperature dependence in perturbative QCD has been calculated to full leading-order for three flavors of massless quarks

 ηs=5.12g4ln(2.42/g) (60)

which can be combined with a two-loop renormalization group expression for the running coupling

 1g2(T)=98π2ln(TΛT)+49π2ln(2ln(TΛT)) (61)

with MeV [49, 50]. From this we take only the most basic result: the shear viscosity to entropy density ratio should increase like the square of the logarithm of the temperature. Assuming that this combines continuously with the shear viscosity at the critical temperature, which we also take as a parameter, we add one more parameter for the slope of increase in shear viscosity to entropy density in the liquid phase

 η/s=η/s|Tc+αln(T). (62)

Large values of have a significant effect on observables and will be discussed in a later chapter.

A reasonable estimate for the bulk viscosity is more difficult to generate. The methods used to estimate the shear viscosity or its scaling properties rely on calculations from models either of massless gases or a collection of mass states where the bulk viscosity is nearly zero. However, near the critical region, bulk viscosity can be finite due to the changing of the sigma mass or effectively be non-zero due to loss of chemical equilibrium [44]. The structure and size of such a peak in the bulk viscosity are not well constrained theoretically and the width, height, and location of this peak are introduced as parameters to the model. These parameters are not explored within this work but are good candidates for future exploration.

## Chapter \thechapter TRISH – A Three-Dimensional Israel-Stewart Hydrodynamics Algorithm

In the first two chapters, we developed hydrodynamic equations of motion and justified their use in the description of the quark matter created in heavy ion collisions. Since hydrodynamics requires an equation of state and transport coefficients from the underlying microscopic theory, these were calculated in the immediately previous chapter. With these results in hand, we can now discuss the algorithm to evaluate viscous hydrodynamics for heavy ion collisions can be discussed. This chapter concerns the algorithm itself including the coordinate system and integration scheme. We then discuss the verification of the algorithm by checking conserved quantities and comparing to known results. Finally, we discuss the process by which the calculation is coupled to the gas calculation that will be used at lower temperatures. This section focuses on the key issues of finding the emission surface and populating particle states once the surface is known.

### 6 Description of Hydrodynamic Algorithm

#### 6.1 Coordinates, Variables, and Integration Scheme

As discussed in Chapter 1, the Bjorken ansatz that the initial longitudinal distribution does not depend on rapidity for small rapidity produces a Hubble-like expansion that hydrodynamics maintains. Even for a finite system where longitudinal density gradients spoil the symmetry, this expansion dominates. For this reason, the algorithm is written in the coordinate system of this expansion. To refresh the memory, lab time and longitudinal position are replaced by proper time and spatial rapidity,

 τ2=t2−z2,  tanhη=z/t. (63)

The integration of the partial differential equations of Israel-Stewart hydrodynamics will take place on an Eulerian grid, meaning that the grid points are at fixed locations in the - coordinate system.

Since the shear tensor is an independent variable in the Israel-Stewart framework, the ten components of the stress energy tensor not given by symmetry are all independent. In principle, the equations of motion could be written in terms of these. Instead, we choose a set of variables more closely related to the original formulation. For instance, ideal hydrodynamics is naturally interpreted in terms of the energy density and three components of collective velocity. From these, the stress energy tensor can then be computed via

 Tμν=(ϵ+P)uμuν−gμνP+Πμν, (64)

where encapsulates both shear and bulk viscous corrections, which would be zero for ideal hydrodynamics. A common alternative in ideal hydrodynamics is to integrate for the temporal part of the stress energy tensor (), which has the benefit of simple equations of motion and a long history of numerical methods to solve conservation equations stably and accurately. [51, 52] However, Equation 64 for is not invertible for , so recovering the collective velocity requires a root find. This is an acceptable cost if it can be done once and stored but doubles the memory cost. This makes it attractive for ideal hydrodynamic codes even with large grids but the larger set of variables make the memory footprint a more important consideration.

Furthermore, while more complex integration methods for conservation equations are an integral part of computational fluid dynamics, it is not clear that they are necessary for a strategy with physical viscosity. Instead we choose a simple and quick integration scheme that is easily multithreaded under the assumption that physical viscosity and smooth initial conditions render such complications unnecessary. Even for these conditions, shear viscosity of around half the entropy density often leads to instability as the correction terms get large. Also, the peak in bulk viscosity near the phase transition can produce shock waves which are not treated by our integration method though for small bulk viscosity the integration is stable.

For the viscous part of the evolution, we choose to track a projected version of the shear tensor in the frame of the matter. This complicates the derivation of the equations of motion as the evaluation of derivatives in this frame requires Lorentz transformations. In addition, the usual geometrical methods of computing corrections for the motion of the coordinate system must be made using small boosts, the details of which are presented in the next subsection. The choice of the local tensor is motivated by the more natural comparison with ideal hydrodynamics – shear corrections can be compared to the pressure without kinematic corrections. This is especially important if one seeks to damp large corrections to the shear tensor in the presence of large velocity gradients where corrections should be restricted to remain smaller than the pressure [35].

The shear tensor has only five independent components due to symmetry, tracelessness, and orthogonality to the collective velocity. Therefore, we take a projection that is inspired by spherical harmonics,

 a1=12(~πxx−~πyy); a2=12√3(~πxx+~πyy−2~πzz); (65) a3=~πxy; a4=~πxz;     a5=~πyz; b = (1/3)[~πxx+~πyy+~πzz].

In this space, is the difference in effective pressure between the x- and y-directions, and is the difference between the longitudinal direction and the average of the transverse directions. Therefore, an infinite transverse system undergoing a Bjorken expansion, only is non-zero, and even when including non-trivial transverse expansion, will be the strongest shear term since it quantifies the extent to which the longitudinal expansion proceeds more rapidly than the transverse expansion. The remaining three still correspond to the shear along the surfaces of the fluid element. Since time derivatives in the Israel-Stewart equations always appear as scaled versions of a fluctuation, we use this scaled anisotropy:

 αi=ai/ση;  β=b/σζ. (66)

The fluctuations are functions only of the energy density, as discussed in Equation 33, and can therefore easily be calculated for each fluid cell at the same time as the other transport coefficients. Since the stress-energy tensor has ten independent elements after symmetry constraints in the viscous case, these six elements combined with the energy density and the collective velocity vector are sufficient to describe the stress-energy tensor completely. The equations of motion are then written in terms of time derivatives of these ten variables exclusively, the derivation of which is the thrust of Section 4.1.4.

It is then convenient to rewrite the velocity gradients to which the shear tensor relaxes in terms of the these projections. Simply applying evaluating these linear combinations of velocity gradients yields

 ω1=~∂xnx−~∂yny, ω2=1√3[~∂xnx+~∂ynx−2~∂znz], (67) ω3=~∂xny+~∂ynx, ω4=~∂xnz+~∂znx,    ω5=~∂ynz+~∂zny (68)

so that Equation 24 can be rewritten as

 τπ~∂tαi=−αi−ησηωi. (69)

#### 6.2 Evaluation of Local Derivatives

The design of our program calls for evaluation of local derivatives of the velocity relative to the expansion of the mesh, , and of the local shear tensor – – where tildes indicate quantities in the local frame. We leave discussion of the expansion corrections to the next subsection and focus on corrections for the motion of the fluid frame. Our goal is to begin with the equations of motion in the fluid frame, and boost such that we have equations of motion for our chosen set of variables in terms of derivatives available in the mesh frame.

We will make frequent use of a general form for the boost from a frame that observes a velocity to one that observes velocity :

 Λμν(u→n)=gμν+2nμuν−(uμ+nμ)(uν+nν)1+u⋅n. (70)

which fulfills the requirement that .

We begin by evaluating the local derivative of the mesh frame velocity. This requires boosting both quantities, , where {1,0,0,0} is the local frame velocity and is the mesh velocity. Considering first the boost to the derivatives

 ~∂μ = Λ νμ∂ν, ~∂μ = [g νμ+2nμuν−(uμ+nμ)(uν+nν)1+u⋅n]∂ν, ~∂μ = ∂μ+2nμuν∂ν−uμ+nμ1+γ(uν∂ν+nν∂ν), (71)

where is the Lorentz factor. Equation 71 is neatly divided into local time derivatives and local spatial derivatives

 ~∂t=γ∂t+ui∂i,  ~∂i=∂i−ui1+γ(uν∂ν+∂t), (72)

where the latin index is used to indicate a spatial index. If is small , this yields the usual comoving derivative expressions

 ~∂t=∂t+ui∂i,   ~∂i=∂i−ui∂t. (73)

Likewise, obtaining derivatives of the local velocities just requires a boost

 ~∂αnβ = Λβν~∂α\emph{u}ν, (74) = [gβν+2nβuν−(uβ+nβ)(uν+nν)1+γ]~∂αuν, = ~∂αuβ−uβ+nβ1+γ~∂αγ, (75)

where we have made use of the fact that

 2uμ∂νuμ=∂ν(uμuμ)=∂ν(1)=0 (76)

assuming that the derivative is a total derivative and coordinate system effects (affine connections) are included. A useful conclusion of Equation 76 is that one can convert derivatives of the Lorentz factor into derivatives of the velocities via

 γ∂μγ=−ui∂μui (77)

for any complete derivative . We will tend to leave derivatives of the Lorentz factor in equations and it is understood that they will be converted into derivatives of the velocities using this relation.

Combining Equations 71 and 75, we obtain an expression for local spatial derivatives of local velocities

 ~∂inj = ∂iuj+ui˙uj+uiuk1+γ∂kuj (78) −(ujukγ(1+γ))[∂iuk+ui˙uk+uium1+γ∂muk],

where summations over the latin indices are implied and over-dots indicate time derivatives in the mesh frame.

The expansion rate, , is a Lorentz scalar and should be conserved under Lorentz transformations. We calculate the local time derivative of the Lorentz factor, , to confirm this. Using 75,

 ~∂0n0=~∂0γ−[1+γ1+γ]~∂0γ=0, (79)

which is expected since, to lowest order in velocity, the Lorentz factor is proportional to velocity squared for small velocities. In fact, nothing about Equation 79 changes for spatial derivatives, which are also zero. If one then sums over spatial indices in Equation 78,

 ~∂ini = ∂iui+[γ−uiui1+γ]∂0γ+[−1−uiui1+γ+γ]11+γuj∂jγ ~∂ini = ∂0γ+∂iui (80)

where the last step makes use of the identity . Combining Equations 78 and 80 confirms that the expansion rate is invariant.

The only remaining velocity derivative to compute is the time derivative of the spatial components, which are

 ~∂0ni=˙ui−ui1+γ˙γ+uj∂jui−uiuj1+γ∂jγ. (81)

For convenience in writing the equations of motion in a subsequent subsection, we will condense the spatial derivatives

 ~∂0ni = ˙ui−ui1+γ˙γ+(~∂0ni)′, ~∂inj = Missing or unrecognized delimiter for \right (82)

Calculating derivatives of the local shear tensor differs somewhat from the previous as the frame of the shear tensor is different for each fluid element. Therefore, we introduce a boost from a frame that observes a slightly different collective velocity to the fluid frame, which we denote . Keeping terms only linear in the difference, we obtain

 δΛμν(u−δu