Chiral symmetry breaking and chemical equilibrium in a heavy-ion collisions

We examine the thermalization of an ensemble of the octet of pseudoscalar mesons, in the isospin symmetric limit, whose interactions are constrained through chiral symmetry, unitarity, and measurements. The reaction amplitudes generate all resonances up to masses of about 2 GeV, with twelve input parameters, namely f_pi, three masses, and eight low energy constants (LECs) of chiral perturbation theory. In linear response theory, we find that matter takes an extremely long time to thermalize. These long relaxation times are directly related to the fact that these mesons are pseudo-Goldstone bosons of chiral symmetry breaking. This result indicates that fireballs created with zero baryon number in heavy-ion collisions will drop out of chemical equilibrium once they enter the chiral symmetry broken phase.


I. INTRODUCTION
In heavy-ion collisions, as in all collider experiments, the main observables are the particles in the final state and their momenta, everything else of interest is constructed from these. The simplest of observables are the abundances of hadrons in the final state. In heavy-ion collisions, these yields are well explained by an ideal gas of hadronic resonances at a freeze out temperature of T f o = 158.4 ± 1.4 MeV [1] when the nucleon-nucleon center of mass collision energies, √ S, is larger than 20 GeV. Furthermore, for √ S > 100 GeV, the flavour chemical potentials are small. Many variants of such models have been examined [2], and they are in reasonable agreement with this result.
A coincidence arises from QCD with its approximate chiral symmetry. In the limit of exact chiral symmetry, there would have been a critical point at finite temperature, and in thermodynamic equilibrium. Since the symmetry is approximate, instead there is a broad cross over [3,4] with a peak in the chiral susceptibility at T co = 156.5 ± 1.5 MeV. This has led to the identification of T f o with T co [1]. At present this is the only point of contact between heavy-ion collisions and the broken chiral symmetry of hadronic physics.
Right from the early days of heavy-ion collisions there have been efforts to build a complete dynamical description of the whole history of the fireball in terms of a transport theory. These transport computations assume that the fireball is made up of an interacting system of quarks and gluons in the very initial stages, trace their interactions and usually find that it approaches equilibrium, cools, and turns into an interacting system of hadrons. In this widely accepted view of the fireball, the number and momentum distribution of the final state hadrons is a result of strong interaction dynamics. T f o is the freezeout temperature, i.e., the point at which the expansion rate of the fireball matches the rate of interactions, and chemistry, and eventually momenta, can no longer be kept in thermal equilibrium through interactions.
Codes which implement these approaches, like URQMD [12] and AMPT [11], incorporate a lot of known physics of strong interactions and try to predict the complete course of a collision. The hadronic interactions which are included in these approaches use many two-to-two hadron processes. A large fraction of these are not measured yet, and have to be constrained by various model considerations. We partly follow this approach in the sense that we examine transport theory with the approximation of two-to-two interactions between hadrons. However, unlike those codes, it is not our aim to model the full course of every collision across a large range of √ S. Instead we ask a more limited question; namely, what would be the relaxation time in a hadron gas pushed slightly out of chemical equilibrium and allowed to relax back to equilibrium. The first step is to decide which hadrons should be included. Certainly pions, which are the lightest of hadrons, need to be accounted for. It turns out that the lightest baryon, i.e., the proton, has an equilibrium number density which is almost two orders of magnitude smaller at T f o . Since pion-proton and pion-pion cross sections are roughly comparable, within the accuracy of a few percent in the relaxation rates, one may neglect the baryon contribution in this computation. In this first study we will do this. Inclusion of strangeness requires us to add kaons, the lightest strange particle, to the mixture. Then the SU(3) flavour structure would require us to add the η. So, the model contains the full octet of the light pseudoscalar mesons. These are the pseudo-Goldstone bosons of approximate chiral symmetry in QCD.
Even though it is approximate, chiral symmetry is predictive because it strongly constrains the low-energy interactions of the pseudoscalar mesons which are the pseudo-Goldstone bosons of this symmetry breaking [6][7][8]. Such a theory is known to be very good at predicting many properties of the lightest pseudoscalar mesons including decay constants and reaction cross sections [8][9][10]. In this work we use the unitarized cross sections of [10].
It turns out that they reproduce the full resonance spectrum, at least up to masses of 2 GeV. So, using the pseudoscalar octet with unitarized amplitudes from chiral perturbation theory seems to capture a large part of the mesonic physics that goes into the resonance gas model, while also giving enough information to make a start on transport computations. Among other pleasant aspects of the computation is that it requires a very small number of input parameters, namely m π , m K , m η , f π , and eight low energy constants of the chiral effective theory (EFT) model. The amplitudes are unitary and contain no further UV cutoffs. The eight LECs are obtained from other hadronic observables. We will keep track of the error bands on all of them.
Some limitations of this computation are clear enough. The neglect of baryons is a major approximation. We plan to remove this in a later work. This should extend the reach of our computations to lower values of √ S. Another major limitation is that the hadronic approach is unlikely to be a reasonable way to capture the physics of the chirally symmetric state of QCD; however, that is not our concern in this work. Our main concern is to obtain a treatment of transport theory in the chiral symmetry broken region of QCD using a controllable, and independently testable, hadron EFT. A step towards this is what we present here.

A. Chemical rate equations and relaxation times
The kinetic theory underlying chemical equilibrium and freezeout is well known [17]. In order to set up our notation and the model approximations, we give a brief review here of the passage from the Boltzmann to the chemical rate equations. The Boltzmann equation for a species a in the reactive fluid can be written in the form where ρ a (x, p) is the Lorentz invariant phase space density, and x and p are 4-vectors for the position and momentum. We write the Liouville operator in the form D = p µ ∂ µ , appropriate for rectilinear coordinates in flat space. The Lorentz vector number current is defined as and g a is the phase space multiplicity factor, which counts, in general, the dimensions of both the spin and isospin representations, i.e., g a = (2S a + 1)(2I a + 1). In our particular application, since we include only pseudoscalar mesons, S a = 0. Note that the particle number density n a is one component of this four vector. In the following, we will examine it in the frame in which the heat bath is at rest. Integrating both sides of eq. (1) one may then write The right hand side of this chemical rate equation can be written as where the sum is over all reactions r(i, f ) which include the particle a among the set of initial particles i and with the appropriate set of final particles f . The reverse reactionr(f, i) interchanges the initial and final sets. Since QCD at zero chemical potentials preserves CP symmetry, one can equate |M | 2 r and |M | 2 r . Further, as long as quantum effects like Pauli blocking or Bose enhancement can be neglected, one can set ρ ≪ 1. As a result, 1 ± ρ ≃ 1. using these two approximations one can write There are similar equations for the whole coupled chain of reactions. No assumptions need to be made at this stage about whether to use Boltzmann or quantum distributions for the distribution functions ρ. This is the form developed in [17], for example.
In the specific case that is of interest here, the temperatures could be slightly higher than the pion mass, but certainly less than twice that; m π < T < 2m π . So reactions which produce more particles in the final state than in the initial are rare. Also, since the phase space densities are much less than unity, collisions of three or more particles are extremely rare. As a result, one can restrict a first investigation to reactions involving two particles in the initial state and two in the final state, i.e., 2-to-2 reactions. In this case, further reduction of eq. (5) is possible.
The total cross section for the forward reaction ab → cd is where the Lorentz invariant definition of the flux of particles in the initial state is , where s is the square of the CM energy in the collision of a and b [18]. One can use a similar expression for the cross section, σr of the reverse reaction, cd → ab. One may trade the squared matrix element on the right for the combination on the left in the rate equation, if one wants to.
We will now choose to work in the frame in which the heat bath is at rest. Our model consists of the fluid at rest in this frame, so that the spatial components of n µ a vanish, and the time component is the particle density. Then the rate equation becomes Here we have used the notation v ab = F ab /(4E a E b ); a Lorentz covariant definition of v ab follows from the above definitions. Further, the double angular bracket indicates averaging over the instantaneous densities of the initial particles for a reaction. However, these non-equilibrium distributions are not universal, and a general study is not of much interest. It is more useful to examine this in the linear-response limit as the system reaches close to equilibrium. Accordingly, we set n a = n eq a + δn a , expand in the small deviations from equilibrium, δn a , etc., and retain terms up to linear order in these small parameters. The terms independent of the δns vanishes due to detailed balance. We can then replace averages with respect to the non-equilibrium distributions by averages in equilibrium. These are denoted by single angular brackets below. To linear order the reaction rate equations become The linear system of equations can be represented in the formṅ = −An, where n is a column vector whose elements are each of the deviations of the number densities of interest from their equilibrium values . Each element of the matrix A is a quantity of the form σv n. One sees that dimensionally this is the inverse of a time (we have used natural units throughout). For every conserved quantity, one has a zero eigenvalue of A. Every other eigenvalue of A is the inverse of one of the relaxation times of the system. The longest relaxation time tells us how fast the disturbed system relaxes back to chemical equilibrium. Far away from equilibrium the system can relax either faster or slower. However, as the system approaches equilibrium, the rate of approach to thermal equilibrium cannot be faster than the inverse of the longest relaxation time. We relate these relaxation times to freeze out in the next subsection.

B. Expansion time scale and freeze out
If chemical freeze out occurs somewhat late in the lifetime of the fireball, then one may approximate the fluid flow by a uniform radial flow. It is most convenient in this case to choose radial coordinates in the CM frame of the fireball. We can then write the Liouville operator in these curvilinear coordinates and find that the linearized form of the chemical rate equations become dn dt where the components of n are the deviations from the equilibrium value of each of the densities, and A is the same matrix as before, and t is the elapsed time in the chosen frame. A small formal remark: the square matrix A is not symmetric. Since the vector n is defined to be a column vector, and A acts on it from the left, the right eigenvectors are the ones which specify the normal modes of the system. The eigenvectors need not be orthogonal to each other. This equation now allows a simple dimensional argument for the freezeout time.
Each eigenvalue of A, λ i , is the inverse of the relaxation time, τ i = 1/λ i , of one of the eigenmodes of the fluid when there is no flow. Every conservation law gives an exact zero eigenvalue; these may be disregarded for the analysis of freeze out. If the other τ i ≪ t/3, then the any deviation from equilibrium dies away much faster than the expansion rate, and the fluid may be considered to be in chemical equilibrium as it expands. If one of the τ i > t/3, then the corresponding eigenmode of A will not be able to relax back to equilibrium, and that mode may be considered to have frozen out. The subsequent evolution of the fluid involves the frozen mode(s) as well as modes which may remain in equilibrium. If there is a big hierarchy between the chemical relaxation times, τ i , then it may be possible to use sequential freeze out scenarios, which are common in cosmology [17], and have been proposed in heavy ion collisions [2].
We note that at very large times, t, the formal solution of the above equation may still allow τ i ≪ t. However, the system will have diluted to the point that a hydrodynamic flow is no longer feasible. Even if the τ i are so small that radial flow is not a very good approximation at t ≈ τ i , a dimensional argument shows that hydrodynamic expansion of the fireball would change the factor 3/t to ξ/t where the dimensionless number ξ depends on the details of the flow. With this change the arguments given above would continue to be applicable.

C. Hadron cross sections from chiral perturbation theory
As discussed earlier, in this paper we report an investigation of a model hadron fluid made of the pseudo-Goldstone bosons of chiral symmetry breaking, namely the lowest SU(3) flavour octet of pseudoscalar bosons. The mutual interactions of these mesons are completely constrained by chiral symmetry [7][8][9], and their reaction cross sections have been computed in chiral perturbation theory. We use the unitarized amplitudes which were presented in [10]. These depend on the meson masses, f π , and 8 other LECs (L 1 , L 2 · · · , L 8 ) which appear in the Lagrangian of chiral perturbation theory to order p 4 .
Discussions of the extraction of the LECs from hadron observables can be found in [8,9]. L 1 , L 2 , and L 3 can be determined through ππ scattering in the J = 2 channel and from the weak decays of K. Constraints on L 1 , L 3 , L 4 , and L 6 come from Zweig's rule. L 5 may be determined by the f π to f K ratio. L 4 and L 6 are also obtained from the π scalar and charge radii. L 5 , L 7 and L 8 are constrained by the Gell-Mann-Okubo mass formulae. Currently error bands in individual constants range from 25-40%. In this work we shall exhibit the uncertainties in predictions due to the current range of uncertainties in these constants.
Two implementations of unitarized amplitudes from chiral perturbation theory are discussed in [10]. In one the dimension of the T-matrix is taken to be equal to the number of channels which are open, so that the dimension changes at every mass threshold. This strictly implements unitarity at each energy, but may result in a discontinuity of the amplitude at thresholds where new channels open up. The other is to use a constant dimension for the T-matrix, equal to the highest dimension required until about 2 GeV. This may result in a loss of unitarity in the vicinity of a threshold, but preserved continuity. The difference between the two methods for phase shifts was shown to be minor. In our implementation we chose the first approach, so that unitarity is always strictly implemented. In Figure 1 we compare our numerical implementation with that of [10]. In the figures we have also given the theory errors on the phase shifts induced by the errors in the input low energy constants. These are estimated by taking the maximum and minimum values of the phase shifts as we vary the parameters along the corners of a hypercube whose center is at the central values of the parameters, and whose faces are one sigma away. We note that keeping track of theory errors is important, as one sees in Figure 1. We also note the spread in measured values of scattering amplitudes from several different experiments. However, some of the errors estimated by moving along the corners of the hypercube are quite extreme, and are dominated by a single parameter. The phase shift δ 1,1 is one such. The magnitude of the discontinuity of the error at the threshold for K pair production is probably an artifact. Because of this we investigate a more natural, but more expensive, estimator. We draw samples of the LECs from independent Gaussian distributions for each parameter with a mean equal to the central value and width equal to the quoted error. Two examples of the error bands obtained by this procedure are shown in Figure 2. Note that the effect of a single outlier is no longer as important. We will use this method in the rest of the paper.
Note also that δ 11 and δ 1/2,1 show very clearly the presence of the ρ and K * mesons. A closer investigation shows that the interacting gas of pseudoscalar mesons generates the full SU(3) octet of vector mesons. We also checked that the interacting system generates all scalar mesons with masses less than 2 GeV. There are no tensor mesons in this mass range, so we the interacting system of SU(3) octet of pseudoscalars generates the hadron resonance gas of all SU(3) octets with masses less than 2 GeV. We use these amplitudes to compute the elements of the matrix A defined after eq. (9).
We note that an earlier work along these lines [5] also tried to use chiral perturbation theory. However, that work predated [10], as a result of which the resonance region was then not captured by this fundamental approach. So, parametrizations were used in [5] to describe amplitudes in the resonance region. We are able to avoid this due to advances in QCD, which allows us to proceed entirely by using a state-of-the-art effective field theory (EFT) with a small number of low energy constants, namely chiral perturbation theory. This allows us to use a small number of parameters which are determined by measurements of other hadron properties, and to write unitary amplitudes for the transport theory without any ad-hoc UV cutoff.

III. RESULTS
In this paper we neglect the differences in masses due to SU(2) isospin breaking, and treat all pions as having a mass equal to the mean mass of the isotriplet. We also take all the kaons to have exactly equal mass, that of the mean mass of each isodoublet. The matrix elements then distinguish four classes of mesons: pions, whose number density we denote by n π , kaons (with strangeness of −1) having number density n K , anti-kaons (strangeness of +1) with number density n K , and η with number density n η . The linearised reaction matrix A is then a 4 × 4 matrix.
Collapsing each isospin multiplet into one species means that we cannot examine thermalization times associated with isospin fluctuations. This is an interesting issue which could be addressed by a fundamental approach such as the one we take. However, the main interest in this question is due to isospin fluctuations in the baryon sector, so we defer this question for later.
The computation of the reaction rate matrix A requires n eq (m, T ). We use the expansion of the number density from the Bose distribution where K 2 is the modified Bessel function of order 2 which decays exponentially as a function of its argument. For temperatures in the range up to 160 MeV or so, ten terms of the series are sufficient to get the number density accurate to five places of decimals for the pion. For the remaining mesons the leading term suffices. This corresponds to the classical Boltzmann gas approximation. The unitarized amplitudes require a partial wave expansion of the expressions obtained from chiral perturbation theory, and, up to the order at which the amplitudes are available, it is sufficient to keep terms only up to J = 2. It is most convenient to do this expansion numerically. Our codes used standard LINPACK routines for unitarization and inversion of matrices, and QUADPACK routines for integrals. Since our matrices are extremely small (3 × 3 at most), we also experimented with simpler high accuracy in-line code for eigenvalues and inverses; these improved run times slightly without compromising accuracy. For the hadronic matrix elements, we found unitarity of the S-matrix in all channels up to machine precision for center of mass energies of a little over 2 GeV. Computations of σv used Gauss-Laguerre routines for the integration over momenta, and Gauss-Legendre routines for angular integrations. We checked that numerical rounding and truncation errors lie well below the theoretical uncertainties due to the input hadron parameters. The comparisons shown in Figure 1 are part of the check of our codes. There are two conserved quantities-net strangeness, and total particle number (since we used only 2-to-2 reactions). As a result there are two eigenvalues of A which vanish. It is simple to actually parametrize the number densities so that this is taken account of. We can write n π = n eq π + h π and n η = n eq η + h η . Then one has n K = n eq K − (h π + h η )/2 and n K = n eq K − (h π + h η )/2. Substituting these parametrizations in eq. (9), one can eliminate two of the equations to obtain the reduced set of equations dh dt where h is a two dimensional column vector whose components are h π and h η respectively, and C is a reduced matrix obtained from A by eliminating the equations for n K and n K using the conservation laws. The eigenvalues of C are the inverses of the relaxation times of the system. The two eigenvectors can be specified by the angle, ψ, that they make with the h π axes. The results of our computations are shown in Figure 3. At a temperature of about 150 MeV, the fast mode has a relaxation time of around 10 fm. This increases to about a 100 fm at a temperature of 100 MeV. The slow mode is an order of magnitude slower, with a relaxation time of around 100 fm at 150 MeV. The slow mode is dominantly of pions relaxing towards equilibrium; the eigenvector makes an angle of 10 to 15 degrees with the pion direction. The fast mode is dominated by the relaxation of the η, since it makes an angle of about 30 degrees with the η direction.
Relaxation times of 100 fm may seem unnatural in a hadron gas, where typical time scales are expected to be around 1 fm. However, one may recall that 1/τ ≃ σn eq (m π , T ), where σ is an average over the hadron cross sections and n eq (m π , T ) is given by eq. (10). The fact that the eigenvector corresponding to the slow mode is almost fully aligned with the pion direction shows that this assumption is fairly accurate. In a theory of pions, one would expect σ ≃ m 2 π /f 4 π . In fact, one finds that the dimensionless ratio Π s = τ s (T )n eq (m π , T )m 2 π /f 4 π (here τ s is the relaxation time of the slow mode) is of order unity, and varies from 4 to 8 over the range of temperature shown in Figure 3. This is a check that the very long relaxation time is natural, and that its magnitude is what should be expected a priori from chiral perturbation theory.
With the relaxation time of the fast mode, τ f , one can similarly compute the dimensionless ratio Π f = τ f (T )n eq (m η , T )m 2 π /f 4 π . Again, this turns out to be of order unity, varying from 0.5 to 1 over the same range of temperature. The use of the number density for η in this case is suggested by the large component of the fast eigenmode in the η direction. Both the fast and slow modes turn out to be natural. The unexpected size of the relaxation times is due essentially to the fact that the equilibrium densities of the mesons at these temperatures are small.
The interactions of Goldstone bosons are strongly constrained, and involve derivative interactions which are forced by symmetry. The fact that quark masses are non-zero allow non-derivative couplings. A signal of this is that the mode dominated by the highest mass pseudo-Goldstone boson, the η, has lower relaxation time. We can check this also by dropping the η meson from the computation. Removing it should decrease the net cross sections, and push up the relaxation time. After removing the two conservation laws using the equations for n K and n K , and writing the equation for the pion in terms of h π , its deviation from equilibrium, one can write the linearized chemical rate equation as Where the relaxation time, τ ′ , of the gas without the η is computed in terms of the remaining rates. The result is about an order of magnitude larger than τ s . The movement is in the direction that we expected from the argument based on chiral perturbation theory.

IV. CONCLUSIONS
We examined the chemical relaxation time in a gas of the SU(3) octet of pseudoscalar mesons in the linear response approximation. The reaction amplitudes in the gas were described by the best current results in chiral perturbation theory available to date [10]. These interactions generated the full SU(3) octet of vector mesons, and all the scalars up to 2 GeV in mass. The amplitudes were completely parametrized by three masses, one decay constant, and eight low energy constants of chiral perturbation theory, all of which were known from experimental observables. We found that the slowest relaxation time, τ s , which controls the rate of approach to equilibrium, is about 100 fm. The second relaxation time, τ f , was found to be of the order of 10 fm. We argued that these long relaxation times are natural for a gas of pseudo-Goldstone bosons. This then implies that such a gas cannot equilibrate chemically within the lifetime of a fireball produced in heavy-ion collisions.
However, this leads to an interesting conjecture. While this study has nothing to say about chemical relaxation times in the chiral symmetric phase of a strongly interacting system, one may assume, as is common, that hot QCD matter has short relaxation times. Then it is entirely possible for this system to come to equilibrium in the early history of the fireball. Our results show that as the system passes through the QCD chiral cross over, T co , and the degrees of freedom change to the pseudo-Goldstone boson of the broken chiral symmetry, it is not possible to maintain chemical equilibrium. This gives a straightforward mechanism for T f o to occur very close to T co . Mild mismatches between the two temperatures are entirely possible due to the fact that a cross over is blurry. This explanation is pleasant because the very mechanism that causes the cross over also causes chemical freezeout.
It should be noted that there are escape routes to this argument. The first is to note that this computation does not include baryons, whereas significant fraction of heavy-ion collision events contain at least a few baryons in the fireball. We are currently investigating a mixture of baryons and mesons, and will report on this aspect later. The next route is technical. The formalism we use is based on a reduction of the Boltzmann transport equations which assumes that the phase space density of particles is significantly less than unity (in natural units). If such deviations from the results of a hadron resonance gas model can be observed, then further investigations will be required. Such an eventuality has been suggested, but has not yet been observed in heavy ion collisions.