Maximum entropy kinetic matching conditions for heavy-ion collisions

Coupling hadronic kinetic theory models to fluid dynamics in phenomenological studies of heavy ion collisions requires a prescription for ``particlization''. Existing particlization models are based on implicit or explicit assumptions about the microscopic degrees of freedom that go beyond the information provided by the preceding fluid dynamical history. We propose an alternative prescription which uses only macroscopic information provided by the hydrodynamic output. This method follows directly from the connections between information theory and statistical mechanics.


I. INTRODUCTION
When modeling heavy-ion collision dynamics macroscopically with fluid dynamics, the problem of particlization of the fluid near its decoupling into hadrons is a persistent source of theoretical model bias in the estimation of the material transport properties of the Quark-Gluon Plasma (QGP) liquid [1][2][3]. Fluid dynamics provides only the hydrodynamic moments of the microscopic distributions of hadrons. For a fluid with conserved energy and momentum (but ignoring conserved charges) the stress-tensor T µν describes the energy and momentum fluxes. It is given by the second momentum-moment of the microscopic distribution, where the four-vectors x and p denote the space-time positions and particle momenta, g h is the spin-isospin degeneracy of hadronic species h, and f h (x; p) is the one-particle distribution function of species h. In local equilibrium the distribution function of each species is uniquely specified by the macroscopic inverse temperature β and the four-velocity of the fluid rest frame u µ . It is given by the Jüttner distribution where θ = 1, 0, or −1 for particles obeying Bose-Einstein, Maxwell-Boltzmann, or Fermi-Dirac statistics, respectively. Out of local equilibrium however, there exist infinitely many microscopic distributions of hadron momenta and yields that reproduce the same hydrodynamic moments. Therefore, practitioners of hydrodynamic phenomenology often choose a particular ansatz for the microscopic physics when particlizing fluid cells in a hybrid hydrodynamic model of heavy-ion collisions. Commonly used ansätze include assumptions regarding the momentumdependence of the viscous corrections to the local equilibrium distribution, such as the Grad ('14-moments') approach [4][5][6], and various approximations of an underlying kinetic theory based on the relativistic Boltzmann equation, e.g. the first-order Chapman-Enskog method with a collision term in Relaxation Time Approximation (RTA) [7,8]. Both the Grad and Chapman-Enskog (CE) methods suppose that the microscopic distribution is split into two terms, the local-equilibrium distribution f eq,h and a dissipative correction δf h , and then solve their respective matching conditions to a finite (usually first) order in the viscous correction δf h (x; p).
In the Grad method one supposes that the viscous correction function may only have a quadratic dependence on the momenta, where b µν are coefficients which are fixed by the matching conditions (1) and thus linearly expressed in terms of the dissipative stresses, the bulk viscous pressure Π and the shear stress tensor π µν . Although at first order this does not involve a microscopic equation of motion [6] (e.g. the Boltzmann equation) it does make a somewhat arbitrary assumption about the possible momentum-dependence of the viscous correction. This assumption can, however, be justified by relating the method to a systematic approximation of the Boltzmann equation in moments of δf [9]. The Chapman-Enskog series with the RTA collision term requires, at leading order, that the dissipative correction satisfies where τ R is the microscopic relaxation time. In practice, it is often assumed that the relaxation time τ R = τ R (x) is species and momentum-independent, and this assumption is implied whenever we refer to the Chapman-Enskog RTA method in this manuscript. Together with the matching conditions (1) this yields an ansatz for the viscous correction δf h which is linear in the dissipative stresses, but with different coefficients than Grad. It was realized early that the viscous corrections from the shear and bulk viscous stresses, motivated by various kinetic theories (the pion gas, perturbative QCD, RTA) had large effects on observables such as the elliptic flow [1,2]. Moreover, large bulk corrections in linearized approaches could overwhelm the equilibrium distribution resulting in the unphysical consequence f h = f eq,h + δf h < 0. This prompted the development of resummation methods designed such that the distribution functions are positive definite for all momenta [10,11].
We note that many of these difficulties faced in the modeling of particlization are not unique to the field of heavy-ion collision phenomenology. Rather, this problem manifests whenever a fluid's coupling is not sufficiently strong, and its expansion sufficiently fast, such that the fluid decouples into particle degrees of freedom. Motivating a microscopic distribution given only knowledge of its macroscopic moments is a generic problem in kinetic theory and statistical mechanics. More broadly, motivating a unique probability distribution given only knowledge of its moments is a generic problem in information theory. The solution (to both) problems that provides the least-biased (maximally-entropic) distribution was given by Jaynes [12]. We follow these methods, tailoring them to the more unique concerns of particlization in heavyion collisions. These methods do not invoke a microscopic equation of motion, nor any ad hoc ansätze regarding the momentum dependence of the viscous corrections. The result is an expression for the microscopic distribution f h that depends non-linearly on the viscous stresses Π and π µν , is positive-definite, matches the entire stresstensor and reduces to the linearized Chapman-Enskog in Relaxation Time Approximation form when the viscous stresses are weak.
Throughout this manuscript we will use natural units h = k B = c = 1 and the mostly-minus metric g µν = diag(1, −1, −1, −1). Lorentz four-vector indices will be denoted by Greek letters while spatial three-vector indices are denoted by Latin letters. Contractions of Lorentz indices will sometimes be denoted by (·), e.g. A µ B µ = A · B or p µ π µν p ν = p · π · p.

II. THE MAXIMUM-ENTROPY DISTRIBUTION
The kinetic entropy density four-current s µ (x) of a system of particles is given by where f h (x; p) is the one-particle distribution function, p the momentum four-vector and x the position fourvector. The function φ[f ] depends on the quantumstatistical nature of the particles through the parameter θ defined in Eq. (2) and is defined by As is the case when performing particlization, let us suppose that some macroscopic theory (e.g. viscous hydrodynamics) provides us with the stress-energy tensor T µν : where is the energy density, p eq the equilibrium pressure, u µ the four-velocity vector, Π the bulk viscous pressure and π µν the shear-viscous tensor. The energydensity and flow velocity are the eigenvalue and timelike eigenvector of the stress tensor The equilibrium pressure can be related to the energy density by an equation of state p eq = p eq ( ), although it will not be necessary to do so in the method proposed in this manuscript. The spacelike projector ∆ µν is defined by It is also convenient to define the symmetric and traceless projector Then, the shear-stress tensor π µν is defined by and the total isotropic pressure is defined by These definitions can be written as constraints on moments of the microscopic distributions f h (x; p): The energy matching condition requires Matching the total isotropic pressure P ≡ p eq +Π requires Finally, matching the shear-stress tensor requires (16) Our approach here is to find the microscopic distribution f h (x; p) which maximizes the entropy density functional given only the ten components of T µν in Eq. (1). We will do so using the canonical method of Lagrange multipliers [12], without imposing any microscopic equation of motion. Our approach differs from Refs. [13,14] where, instead of the entropy, the entropy production rate was extremized. Computing that rate requires a microscopic approach -in particular, one must specify the collision term. This incorporates additional information that we pretend not to possess -in our work, we assume that all that is known is the energy-momentum tensor (1) resulting from the preceding hydrodynamic evolution of the fluid. In the absence of shear and bulk viscous stresses, our approach recovers the local equilibrium distribution (2) -the present work generalizes it to a maximum-entropy distribution for systems with non-zero viscous stresses.
After introducing Lagrange multipliers with the appropriate tensorial structure, the entropy density fourcurrent can be written The Lagrange multipliers Λ = Λ(x), λ µ = λ µ (x) and γ µ αβ = γ µ αβ (x) are all functions of spacetime x, although we do not explicitly write it for brevity of notation. In the last line we observed that π αβ = ∆ αβ ρσ π ρσ ≡ π αβ , where ∆ αβ ρσ is defined in Eq. (11), and used this to simplify the tensor structure of the Lagrange multiplier γ µ αβ . We seek the distribution f h which maximizes the entropy density in the local rest frame u · s: It follows that and, after simplification and exponentiation, yields the main result: We refer to Eq. (20) as the maximum-entropy (ME) distribution. Here p α = ∆ αβ p β denotes the spatial components of p µ in the local rest frame (LRF), defined by u µ LRF = (1, 0). In that frame the maximum-entropy distribution is given by where p 0 = m 2 +p 2 , γ 0 ij is traceless and symmetric in the spatial indices (ij), and we have defined λ Π ≡ u · λ.
This equation can be rewritten in a form which bears resemblance to previous particlization ansätze which have been studied in the past [10], where the linear-transformation operator Λ ij (which acts on the spatial momenta in the LRF) is defined by This maximum-entropy distribution, in particular the tensor structure of the transformation Λ ij , indeed bears striking resemblance to the so-called "modified equilibrium" distributions [10,15]; however, we find the coefficients in Λ ij are different. Eq. (22) also shares some structural similarities with "anisotropic equilibrium distribution" functions [16][17][18][19]; again, a closer comparison reveals differences. We show in section III that if one works to linear order in the dissipative stresses π µν and Π, the maximum-entropy prescription matches exactly the Chapman-Enskog RTA method. This feature is also shared by particular modified equilibrium approaches [10,15], meaning the relation between the maximum-entropy, modified equilibrium and Chapman-Enskog RTA approaches is exact at linear order in the dissipative stresses. However, once second-order and higher terms have been included, this equivalence is broken and all three prescriptions differ.
We note that for Maxwell-Boltzmann particles (θ = 0) the non-equilibrium entropy density of our system admits a thermodynamic expression in terms of the hydrodynamic fields and their conjugate variables. In this case, the entropy density in the local rest frame s is given by where n is the particle density and we've defined the pressure tensor P ij by which contains both the equilibrium pressure p eq δ ij and viscous corrections. A more general thermodynamic relation which holds for Bose-Einstein and Fermi-Dirac statistics as well can be defined by introducing a generating function Z. The expression is given by where the generating function Z is defined such that its derivatives generate the hydrodynamic fields, Our distribution function is expressed in terms of seven unknown Lagrange multipliers Λ and Λ ij (or, covariantly, Λ, u · λ, and γ µν ≡ u α γ α µν ) which must be chosen to match the energy density, total isotropic pressure and shear-stress tensor, respectively. To solve for these coefficients, we write down the seven required matching conditions (14)- (16) in terms of LRF momenta and components as follows: Note that, since we only introduced a single Lagrange multiplier for each constraint, i.e. we maximized the LRF entropy density only subject to the information given to us directly through the energy momentum tensor, the maximum entropy distribution (20) automatically distributes the shear and bulk viscous flows "democratically" [20] across the hadron species h.

III. LINEARIZING THE MAXIMUM ENTROPY DISTRIBUTION
In this section we solve for the Lagrange multipliers in the limit of small dissipative flows. We calculate the maximum entropy distribution self-consistently to leading order in the viscous stresses Π and π µν , and find that it reproduces exactly the Chapman-Enskog distribution in the Relaxation-Time Approximation. We begin by noting that for vanishing viscous stresses (ideal fluids) the constraint (29) is solved by setting Λ ij ≡ 0 and Λ = β, the equilibrium inverse-temperature. The leading order expressions can therefore be obtained by expanding the maximum entropy distribution (20) to linear order in Λ ij and the difference Λ − β. For simplicity we will assume Maxwell-Boltzmann statistics (θ = 0) throughout this section. We also introduce the compact notation where p 0 = p 2 + m 2 h . (Note that the shorthand p includes a sum over species h.) Let us consider Eq. (20) and for notational convenience define the traceless and purely spatial (in the LRF) rank two-tensor γ αβ (x) = γ αβ (x) ≡ u µ (x)γ µ αβ (x). The maximum-entropy distribution (20) is then expressed as This distribution depends on the hadron species h only through the mass dependence of the on-shell energy. The local-equilibrium distribution is given by Let us now consider the relation between Λ, λ Π , γ µν and β provided by the energy matching condition Eq. (29) discussed at the end of the preceding section and expand it around λ Π = γ µν = 0: Here Λ(β, 0, 0) = β, the coefficients c λ and c µν are in general functions of β, and the ellipses denote neglected terms of second or higher order in λ Π and γ µν . To linear order in the Lagrange multipliers the distribution function f h can thus be written The matching conditions for the shear and bulk stresses can be compactly written as where δf h is defined in Eq. (3). Using the linearized form (36) and introducing the shorthand To calculate the left hand side we use the following tensor decomposition for n = 4: p p µ1 p µ2 · · · p µn (u · p) r f eq,h = I r (n,0) u µ1 u µ2 · · · u µn + I r (n,1) ∆ µ1µ2 u µ3 · · · u µn + ∆ µ1µ3 u µ2 u µ4 · · · u µn + permutations + · · · where the moments I r (n,q) are defined as Note that the authors of Ref. [8] used a similar definition, albeit for a single hadron species. Labeling the moments used in [8] as I r,h (n,q) for a given species h, the two definitions are simply related by I r (n,q) = h I r,h (n,q) . After tensor decomposition we find In the last line we used that γ µν is symmetric, traceless and transverse to the flow velocity u µ . Substituting this back into Eq. (38) yields Using the mutual orthogonality of the tensors u µ u ν , ∆ µν and π µν we find the following three relations: With the help of Eqs. (8,11,12,22) in Ref. [8] for a single hadron species, remembering that here I r (n,q) = h I r,h (n,q) as well as β π = h β π,h (where β π,h ≡ βI 1,h (4,2) was defined in [8]), we find In the last equation we used dp eq = h dp eq,h and d = h d h , as well as dp eq,h = I 0,h (3,1) . Putting everything together we find the following relations between the Lagrange multipliers, coefficients and dissipative stresses: where, similar to the corresponding definition in [8], Using these results in Eq. (36) yields for the linearized maximum-entropy viscous correction δf h the expression For a single hadron species this matches exactly with Eq. (27) in Ref. [8] which was derived by solving the first-order Chapman-Enskog correction with a relaxation-time approximation collision kernel.

IV. MATCHING WITHOUT SHEAR STRESS
We now consider the case in which we match an energymomentum tensor with vanishing shear-stress, π µν = 0.
This implies that the associated Lagrange multiplier γ µν = 0 and that the distribution (20) is isotropic in the LRF. Using LRF momenta (i.e. p 0 = u · p) to evaluate the matching integrals (29)-(30) we obtain FIG. 1. The particle spectra in the local fluid rest frame for the local equilibrium (solid), maximum-entropy (ME, dashed) and linear Chapman-Enskog RTA (CE, dotted) distributions, for pions (blue) and protons (red). The bottom panel displays the ratio of these spectra to the local equilibrium spectra.
Examining the expression for the distribution function where Λ > 0, we note that existence of a solution requires λ Π ≥ −Λ. Since large values of |λ Π | signal large bulk viscous stresses, we restrict our solution to the case |λ Π | ≤ |Λ|. 1 To solve Eqs. (49,50) numerically we use a realistic hadron resonance gas (HRG) which includes all resonances that can be propagated in the UrQMD [21,22] hadronic afterburner. We solve these coupled equations for Λ and λ Π as follows: First, we choose a regular grid of values for both Λ and λ Π . Then, for each pair of values (Λ i , λ Π,j ) we evaluate the integrals in Eqs. (49) and (50) by numerical quadrature. This yields a grid of values (Λ i , λ Π,j ) and P (Λ i , λ Π,j ). These grids are then interpolated with splines to obtain smooth approximationŝ e(Λ, λ Π ) andP (Λ, λ Π ). Finally, given known values of energy density and isotropic pressure, Λ and λ Π can be found using a two-dimensional root finding routine.
We now compare the maximum-entropy distribution with the linear Chapman-Enskog RTA distribution and the local-equilibrium distribution, all in the local rest frame. The linear Chapman-Enskog RTA viscous correction is given by [8,11,15,23] δf where the coefficients F and β Π are given by thermodynamic integrals of the hadron gas: where the moments J (n,q) are defined by (54) For our comparison we consider a hadron resonance gas at temperature T = 0.15 GeV. We assume a moderately large, negative bulk pressure with magnitude of one-third of the hadron resonance gas equilibrium pressure at this temperature: Π = −p eq /3. For a bulk pressure of this magnitude the Chapman-Enskog RTA linear bulk correction becomes larger than the ideal part, |δf h,CE | > ∼ f h,eq , already at moderate values of momentum |p| ∼ 1 GeV. In practice, when doing particlization in simulations of heavy ion collisions, the viscous correction must be regulated by hand. This is required to maintain positivity of the distribution function, which is interpreted as a probability distribution from which particles and their momenta are sampled. A typical procedure is to replace which we will call the 'regulated Chapman-Enskog viscous correction'. We note that, in practice, this regulation breaks the exact matching of the dissipative part of the stress-tensor: Restoring the exact matching condition would require recalculating the coefficients F and β Π using the regulated Chapman-Enskog RTA correction form in the integrands. In Fig. 1 we compare the maximum-entropy distribution with the regulated CE RTA distribution and the local equilibrium distribution. The deviations from the equilibrium distribution are large. For pions, the viscous corrections are negative, except at very low momenta p < 50 MeV. The CE RTA and ME distributions agree well at low momenta but disagree dramatically at p > 1 GeV where the (unregulated) CE RTA distribution goes negative. For the heavier protons, the bulk viscous corrections are much larger and switch sign at intermediate momenta (p ∼ 0.7−0.9 GeV), being positive at lower and negative at larger momenta. Significant discrepancies between the two viscous distributions are observed for protons over the entire momentum range; at p = 0, the ME proton distribution is about 20% larger than the CE RTA distribution and more than a factor 2 larger than the equilibrium distribution.
In Fig. 2 we plot the local rest frame particle densities and mean momentum magnitudes, both normalized by their equilibrium values, for pions and protons as functions of the bulk inverse Reynolds number −Π/p eq . We see that for protons the maximum-entropy distribution yields smaller bulk viscous corrections to the particle density than the regulated linear Chapman-Enskog RTA result. For pions the bulk viscous corrections to the particle densities have the opposite sign and are much larger, but the predictions of the ME and regulated CE RTA distributions agree well with each other, even for large values of the bulk viscous pressure. The same is not true for the pion and proton mean momenta which, for large bulk inverse Reynolds numbers, differ significantly, in opposite directions, between the maximum-entropy and regulated Chapman-Enskog viscous distributions. The differences in the proton yields and pion and proton mean momenta between these two different ansätze for the bulk viscous distribution functions are large enough to have the potential of significantly affecting the bulk viscosity inferred from model-to-data comparisons.
As a final note we observe that the maximum-entropy distribution can naturally handle very large bulk inverse Reynolds numbers. Physically, the bulk inverse Reynolds number may grow large near the pseudo-critical temperature T c ∼ 150 MeV where the quark-gluon plasma turns into hadrons [24][25][26]. Since the maximum-entropy method is not based on a near-equilibrium expansion, it does not require the dissipative stresses to be small for self-consistency.

A. The relationship between the shear stress and its Lagrange multipliers
In this section, we again consider the simpler case of Maxwell-Boltzmann particles (θ = 0). We return to the matching condition (16) for the shear stress tensor and write it in the local rest frame: (57) Latin tensor indices run over the spatial directions in the LRF.
The matching condition (57) establishes a highly nonlinear relationship between the given shear stress in the LRF, π ij , and the associated symmetric and traceless 3×3 matrix of Lagrange multipliers, γ ij . We now proceed to show that these two matrices share a common set of eigenvectors. This will be seen to decisively simplify the task of determining the Lagrange multipliers γ ij from the shear stress π ij .
Let us Taylor expand the second exponential in (57), and consider truncating this series at some finite order. Truncating at n = 0 yields π ij (0) = 0. At truncation order n = 1 we get where C 1,0 (Λ, λ Π ) is a scalar function. In C r,s , the first index r denotes the total number of tensors γ following the coefficient while the second index s denotes how many of these γ tensors have their indices mutually contracted to form scalars. (This will become clearer below.) At n = 2 we encounter the integral ∆ ij kl γ ab γ cd where (γ 2 ) kl = γ k c γ cl . (Note that the term C 2,1 γ kl γ a a vanishes because γ is traceless.) The integral over the n = 3 term in the series yields where we introduced γ n ≡ tr(γ n ), and the n = 4 term similarly integrates to ∆ ij kl C 4,0 (γ 4 ) kl + C 4,2 (γ 2 ) kl γ 2 + C 4,3 γ kl γ 3 . (62) It is clear that the n th order introduces one new term while all other terms have the same tensor structure as lower order terms. For example, when we truncate Eq. (58) at order n = 4 we obtain with the coefficients c (4) The scalar coefficients c i are functions of Λ, λ Π , and scalar contractions of γ. Summed to all orders the series for π ij thus becomes where the coefficientsc n include contributions c (i) n from all orders 1 ≤ i ≤ ∞. Note that although γ ij is traceless, products of γ ij 's are not. The operator ∆ ij kl projects out the trace part of such products. For example, Equation (65) can be written in matrix notation as where I is the 3 × 3 identity matrix and Since any power of a matrix commutes with itself, [γ, γ n ] = 0, it follows that [γ, π] = 0. Therefore, the shear stress tensor π and its associated tensor of Lagrange multipliers γ are simultaneously diagonalisable. At any spacetime point x, the hydrodynamic energymomentum tensor T µν (x) provides us (after transformation to the LRF at x) with the full matrix π. Finding the eigenvalues and eigenvectors of the given π is straightforward. As a real, symmetric and traceless matrix, π has two independent real eigenvalues, associated with three orthonormal eigenvectors. [It is easy to show that orthogonality and normalization allow to characterize the three eigenvectors by three real parameters (Euler angles). Together with the two independent eigenvalues we thus recover the five independent shear stress degrees of freedom.] The matrix C of orthonormal eigenvectors can be used to rotate π into diagonal form: Here π D ≡ diag(π 1 , π 2 , −π 1 −π 2 ) is the diagonalized shear stress, with eigenvalues π 1 , π 2 , and π 3 = −(π 1 +π 2 ). From the arguments above it follows that γ is diagonalized by the same transformation: (70) This implies that three of the five independent Lagrange multipliers can be determined easily from the eigenvectors of the shear stress. Since the rotation matrix C is known from the diagonalization (69) of the given shear stress π, the full Lagrange multiplier matrix γ is easily computed from (70) once the two independent diagonal elements of γ D have been determined. The only nonlinear part of the problem is the computation of the two independent eigenvalues of γ from those of π. These statements hold irrespective of the value of the bulk viscous pressure. They greatly simplify the numerical task of finding the Lagrange multiplier matrix γ that matches the given shear stress π. 2

B. Solution for a massless gas
In this Section we illustrate the determination of the Lagrange multipliers γ ij in the maximum-entropy distribution (20,21) from π ij for the case of a massless Maxwell-Boltzmann gas with nonzero shear stress but vanishing bulk viscous pressure. By keeping the discussion initially general we show that the future generalization of this solution to a general gas mixture of massive hadronic resonances characterized by non-zero values for both shear and bulk viscous stresses will be straightforward.
Let us return to the generating function Z from Eq. (27), now expressed covariantly and specifically for Maxwell-Boltzmann particles: (71) 2 The idea to simplify the solution of the matching conditions by computing the hydrodynamic moments of the anisotropic distri-bution function in a frame which diagonalizes the shear tensor was also exploited in Ref. [19].
According to Eq. (28), the energy density, total isotropic pressure and shear stress from Eqs. (14)- (16) are then given by the following derivatives with respect to the Lagrange multipliers: Note that in the generating functional (71) we do not impose tracelessness and transversality on γ αβ , i.e. when taking the derivatives (72) we consider, in particular, all three eigenvalues γ 1 , γ 2 , and γ 3 as independent. In the last equation (72), the correct symmetries of π µν are ensured by the spatial and transverse projection implied by the angular parentheses in the factor p α p β in Eq. (71). As before, we will work out these expressions in the local rest frame. However, to simplify the last factor in (71) we will rotate the LRF integration momentum variables with the matrix C that diagonalizes γ. The integration measure p and the first three factors under the integral (71) are invariant under this rotation. Writing q = C p as well as q ≡ | q|, q 0 = q 2 +m 2 , and computing (the tracelessness condition γ 3 = −(γ 1 +γ 2 ) will only be implemented at the end), the generating function (71) takes the form Here (θ q , φ q ) are the polar and azimuthal angles of q, with standard integration measure dΩ q = sin θ q dθ q dφ q . From Eq. (74) the generating function Z can be calculated numerically. Taking its derivatives (72) with respect to the four independent Lagrange multipliers on which it depends and solving the resulting matching conditions iteratively may be a difficult task whose full solution will be left to future work. In the following two subsections we work out the generating function semi-analytically for the simpler case of a massless Boltzmann gas (where m=0 and q 0 = q = | q|) and then evaluate the matching conditions in this simplified setting.

Alternate method
The following approach uses a series expansion that allows to perform the angular integrals even in the general case of a gas mixture of massive hadrons. Writing the exponential under the angular integral in Eq. (74) as where q 1 = q sin θ q cos φ q , q 2 = q sin θ q sin φ q , and q 3 = q cos θ q , the angular integrals over each term can be done: A(n 1 , n 2 , n 3 ) = 2 Γ(n 1 + 1 2 )Γ(n 2 + 1 2 )Γ(n 3 + 1 2 ) Γ(n s + 3 2 ) , where we introduced n s ≡ n 1 +n 2 +n 3 for brevity. This yields the following form for the generating function: A(n 1 , n 2 , n 3 ) This expression is still valid for a massive gas mixture with classical Boltzmann statistics and can thus form the basis for future generalizations of what we derive below. Please note that we have not yet used the zero-trace condition γ 1 +γ 2 +γ 3 = 0; as before, we only implement it at the end. We now restrict our attention again to the simpler case of massless particles (q 0 = q) for which the momentum integral is easily performed: HereΛ = Λ + λ Π − 1 3 (γ 1 +γ 2 +γ 3 ). By keeping a sufficient number of terms in the series given above, we have checked that both Eq. (82) and Eq. (78) yield identical results for Z for a given set of Lagrange parameters. The derivative with respect to the eigenvalue γ 1 (γ 2 ) of the Lagrange multipler tensor γ yields the eigenvalue π 1 (π 2 ) of the shear stress tensor: The energy density and total isotropic pressure are given by All of these expressions are to be evaluated at γ 3 = −(γ 1 +γ 2 ) andΛ = Λ + λ Π . The series in the last two equations are seen to be related by P = /3 which (with the equation of state for a massless gas p eq ( ) = /3) yields Π = 0, as it should: in a massless gas the bulk viscous pressure vanishes. Correspondingly, λ Π = 0 for this system.
The algorithm and root solver were tested in a blind test where one of the authors selected a value of Λ and a matrix γ of shear stress Lagrange multipliers in the LRF, without any symmetry restrictions, used Eqs. (29) and (31) (with λ Π = 0) to generate the corresponding energy density and the full LRF shear stress tensor π ij for a massless Boltzmann gas, and handed these to another author who then diagonalized π ij , used Eqs. (83)-(85) to reconstruct (Λ, γ 1 , γ 2 ) and, finally, the complete Lagrange multiplier tensor γ from Eq. (70). The reconstructed Lagrange multipliers agreed with the originally selected ones to a precision that can be systematically improved by truncating the numerical series at higher orders. 5 In Figs. 3 and 4 we plot the momentum distribution in the local rest frame for two simple cases. In both cases, we take an energy density which corresponds (for a massless Boltzmann gas) to an equilibrium temperature of T sw = 0.15 GeV.
In the first case, we assume the shear stress is isotropic in (x, z) but anisotropic in (x, y), taking π xx = π zz = p eq /5 (where p eq is the equilibrium pressure) and offdiagonal components zero. It follows that π yy = − 2 5 p eq . Fig. 3 shows the resulting variation of the maximumentropy distribution as a function of the azimuthal angle φ p for particles with momenta of average thermal magnitude p = 3T sw . Also shown is the linear Chapman-Enskog RTA distribution, with viscous correction given by Eq. (48).
FIG. 4. The maximum-entropy (ME, blue dashed) and linear Chapman-Enskog RTA (CE, red dotted) distributions in the local rest frame as functions of polar momentum angle θp (top) and momentum magnitude (bottom) for a diagonal shear stress with (x, y) isotropy, taking π xx = π yy = peq/5. Both are normalized by the equilibrium distribution.
In the second case, we take a shear stress tensor which is again diagonal but now isotropic in (x, y), given by π xx = π yy = p eq /5. It follows that π zz = −2p eq /5. In Fig. 4 the maximum-entropy and linear Chapman-Enskog RTA distributions are plotted for particles with average thermal momentum p = 3T as a function of polar angle (left panel), and for fixed direction θ p = π/2, φ p = 0 as functions of the momentum magnitude (right panel).
FIG. 5. The normalized Lagrange multiplier γ1/β (where β is the inverse equilibrium temperature) as a function of the first shear stress eigenvalue π1, for two choices of its second eigenvalue π2, π2 = 0 (blue dashed) and π2/peq = 0.3 (green dashed), respectively. Also shown for comparison is the corresponding CE RTA coefficient (orange dotted) which is independent of π2.
Finally, in Fig. 5 we plot the evolution of the Lagrange multiplier γ 1 (the first eigenvalue of γ) with the associated shear stress eigenvalue π 1 /p eq (which is related to the inverse shear Reynolds number), for two choices of the second shear stress eigenvalue, π 2 /p eq = 0 and 0.3, respectively. It is compared with the corresponding CE RTA coefficient which is linear in π 1 and independent of π 2 . In both cases the energy density is fixed by the equilibrium energy density of the massless Maxwell-Boltzmann gas at a temperature T = β −1 = 0.15 GeV. For π 2 = 0, the two coefficients agree very well, even for large inverse Reynolds numbers. However, agreement in the coefficients should not be interpreted as agreement in the predicted distributions: their functional forms are different (exponential momentum dependence in the ME distribution, polynomial dependence in the linearized RTA CE distribution). The green dashed line in Fig. 5 illustrates that for non-zero π 2 the coefficients γ 1 /β and π 1 /(2β π ) differ and do not agree with each other even for very small π 1 . This is a direct manifestation of the non-linear coupling between the eigenvalues γ 1 and γ 2 in the ME matching of the shear-stress.
Although the method for finding the Lagrange multipliers (Λ, γ µν ) from a given energy-momentum tensor with shear stress were here demonstrated numerically only for a massless Maxwell-Boltzmann gas, the framework for doing so for a general gas mixture of massive hadron resonances with Boltzmann statistics has been provided in this work, and its generalization to account in the generating function Z for Bose-Einstein or Fermi-Dirac statistics should be straightforward. However, a more efficient numerical routine for evaluating the momentum integrals in the massive particle case needs to be developed.

VI. CONCLUSIONS AND OUTLOOK
We have worked out the maximum-entropy distribution function as an alternative prescription for particlizing a fluid in a heavy-ion collision. For a general gas mixture of massive hadron resonances, we were able to solve numerically for the distribution function in the case that there was a non-zero bulk viscous pressure while the shear stress vanished. By comparing with the linear Chapman-Enskog RTA prescription we found that the maximum-entropy method yields significantly different particle momentum distributions and yields which can have non-negligible consequences for the theoretical interpretation of experimental data. For a gas of massless Maxwell-Boltzmann particles we demonstrated an algorithm for finding the maximum-entropy distribution when particlizing a fluid with vanishing bulk viscous pressure but non-zero shear stress. A full numerical solution of the maximum-entropy distribution for a massive hadron resonance gas in which both bulk and shear viscous stresses are nonzero is outstanding but of high value for phenomenological modeling of heavy-ion collisions.
Although in the present work we have not included any conserved charges such as net baryon number and strangeness, the generalization of the maximum-entropy prescription to include related dissipative effects (such as non-vanishing baryon and strangeness diffusion currents) should be straightforward. In general, this method allows to match the distribution function at particlization to any macroscopic quantity of which we have prior knowledge on the particlization surface.
Modern phenomenological studies of experimental heavy-ion collision data aim at reconstructing from the data, with quantified uncertainties, key parameters characterizing the evolving hot and dense medium (see, e.g., Refs. [3,[28][29][30] for very recent examples of this type of approach). This is done within a Bayesian statistical framework in which the inferred probability distribution for the model parameters of interest (the "posterior") is obtained as the product of a "prior" distribution for the parameters (accounting for any prior knowledge that we might possess before performing the model-data comparison) and a "likelihood" which accounts for how well, for a given parameter choice, the model predictions agree with the measurements.
An important consideration in Bayesian inference is to avoid introducing uncontrolled physics models in the likelihood that bias the parameter estimates. If assumptions made about the microscopic physics are not well-justified, the resulting model parameters won't be either. The maximum-entropy distribution introduced in this work provides a functional form for the unknown distribution of particles that implements all of, and only the information given to us by the hydrodynamic theory describing the dynamical evolution preceding the particlization process. In this sense it is the least biased choice that can be made in the absence of a trustworthy microscopic theory of the hadron gas close to the pseudo-critical temperature. Any other choice introduces additional information ("theoretical prejudice") into the particlization process that, as far as we know, cannot be compellingly justified theoretically.

ACKNOWLEDGMENTS
We thank Michael McNelis for very useful discussions regarding the comparisons between the maximum entropy and linear Chapman-Enskog RTA prescriptions. This work was supported by the National Science Foundation (NSF) within the framework of the JETSCAPE Collaboration under Award No. ACI-1550223. Additional partial support by the U.S. Department of Energy (DOE), Office of Science, Office for Nuclear Physics under Award No. DE-SC0004286 and within the framework of the BEST and JET Collaborations is also acknowledged.