Resonance Recombination Model and Quark Distribution Functions in the QuarkGluon Plasma
Abstract
We investigate the consequences of spacemomentum correlations in quark phasespace distributions for coalescence processes at the hadronization transition. Thus far it has been proved difficult to reconcile such correlations with the empirically observed constituent quark number scaling (CQNS) at the Relativistic HeavyIon Collider (RHIC). To address this problem we combine our earlier developed quark recombination model with quark phasespace distributions computed from relativistic Langevin simulations in an expanding QuarkGluon Plasma (QGP). Hadronization is based on resonance formation within a Boltzmann equation which recovers thermal equilibrium and obeys energy conservation in the quarkcoalescence process, while the fireball background is adjusted to hydrodynamic simulations of semicentral AuAu collisions at RHIC. To facilitate the applicability of the Langevin process, we focus on strange and charm quarks. Their interactions in the QGP are modeled using leadingorder perturbative QCD augmented by effective Lagrangians with resonances which smoothly merge into hadronic states formed at . The interaction strength is adjusted to reproduce the empirical saturation value for the quarkelliptic flow, . The resulting and elliptic flow recover CQNS over a large range in transverse momentum () within a few percent. As a function of transverse kinetic energy, both the quark spectra from the Langevin simulations and the meson spectra generated via resonance recombination recover CQNS from zero to at least .
I Introduction
The mechanism of hadronization, i.e., the conversion of quarks and gluons produced in hadronic or electromagnetic reactions into colorless hadrons, is a nonperturbative problem that is presently not calculable within the theory of strong interactions, Quantum Chromodynamics (QCD). For partons produced at large transverse momentum, , the factorization theorem of QCD allows to treat the hadronization process via a socalled fragmentation function, which is universal and can, in principle, be determined empirically. At low momentum this scheme is no longer applicable and other hadronization mechanisms become relevant. The quark coalescence model (QCM) has provided a phenomenologically successful framework to understand several nonperturbative features of hadron production in hadronic collisions. In elementary ( and ) reactions, flavor asymmetries in kaon and charmed hadron spectra have been associated with the recombination of produced strange and charm quarks with valence (and sea) quarks in target and projectile (1); (2); (3). In heavyion reactions at the Super Proton Synchrotron (SPS) (4) and the Relativistic HeavyIon Collider (RHIC), recombination of quarks from a thermalized QuarkGluon Plasma (QGP) (5); (6); (7); (8); (9); (10); (11); (12); (13); (14) gives a simple and intuitive explanation of several unexpected features in the observed hadron spectra, most notably the large baryontomeson ratio and the rather universal constituent quark number scaling (CQNS) of the elliptic flow coefficient, (: number of valence quarks in hadron , : transverse momentum of ). This scaling relation implies the momenta of the coalescing quarks to be collinear, and its implementation in QCMs usually restricts their applicability to sufficiently large momenta so that the associated nonconservation of energy in the hadron formation process is small. Since at high parton fragmentation is expected to take over, the typical range of applicability of QCMs is at intermediate momenta, . In our recent work (12), we have suggested a reinterpretation of quark coalescence in terms of hadronic resonance formation, implemented via scattering into a Boltzmann equation (: meson). Energy conservation is obeyed by utilizing hadronic reaction rates, along with detailed balance, based on pertinent spectral functions. In addition, we have shown that this approach correctly recovers the thermal equilibrium limit, which enabled a more controlled extension of the coalescence mechanism to low and to make contact with the phenomenologically successful hydrodynamic description of bulk matter at RHIC.
Another aspect that has evaded a satisfactory explanation in QCMs at RHIC is the question of spacemomentum correlations in the underlying (thermal) quark distribution functions (see, e.g., Ref. (15) for a recent critical review). In hydrodynamic models the elliptic flow of produced hadrons is a collective effect that implies a definite correlation between the particle’s momentum and its spatial position in the fireball, i.e., a (locally thermalized) fluid cell moving into a specific direction preferentially emits hadrons in that same direction. Such a correlation is neglected within the socalled “factorized” implementation of the parton which does not carry any spatial dependence (and is therefore identical regardless of the quark’s position inside the fireball). While this approximation straightforwardly recovers the empirical constituentquark number scaling (CQNS) of the hadron elliptic flow, it is at variance with the hydrodynamic description of as a collective expansion effect.
Previous attempts to incorporate spacemomentum correlations into QCMs have found the empirically observed CQNS to be rather fragile (16); (17); (18). Part of the problem is the construction of a realistic transition from the thermal to the kinetic regime of the underlying parton phasespace distribution functions, as characterized by the “saturation” (leveling off) of the empirical parton at about . In Ref. (16) several elliptic “deformations” of a thermal blastwave parameterization have been considered, motivated by different plausible realizations of . While some features could be ruled out being incompatible with the empirical CQNS, other assumptions did not spoil the latter. In Ref. (18) a reduction of the boost velocity at higher momenta was introduced, entailing a violation of CQNS at the 20% level. It therefore seems that purely phenomenological prescriptions of spacemomentum correlations and associated did not arrive at a conclusive interpretation of the key features underlying the quark distribution functions. In Ref. (17), based on numerical transport simulations, it was even argued that rather delicate cancellations must be at work to obtain CQNS for the thermal components, thus raising doubts on the robustness of the coalescence approach. In view of the broad empirical applicability of CQNS across different centralities, system sizes and collision energies (19); (20), such an interpretation would be difficult to reconcile with experiment.
In the present paper we adopt a microscopic approach to compute quark distributions in fourdimensional phase space (transverse position and momentum) by employing relativistic Langevin simulations for strange and charm quarks within an expanding thermal QGP background. In a strict sense, the underlying FokkerPlanck equation is applicable for a diffusive treatment of heavy and/or highmomentum quarks, i.e., in a regime where the momentum transfers from the heat bath are small. Our simulations for lowmomentum () strange quarks are thus at the boundary of applicability of a FokkerPlanck treatment and may be considered as extrapolations thereof. The Langevin approach has the attractive feature that it naturally encodes the transition from a thermal to a kinetic regime. In particular, when simulating heavyquark (HQ) motion in an expanding QGP fireball for noncentral collisions, this transition reflects itself in a saturation of the elliptic flow (21); (22), a key ingredient to CQNS in lighthadron spectra observed at RHIC. It turns out that Langevin simulations preserve the saturation feature when applied to strange quarks (with thermal masses of ). We will therefore investigate whether the resulting quark distribution functions, evolved to the hadronization transition and injected into our resonance recombination approach, allow for a better (microscopic) understanding of spacemomentum correlations in the coalescence process. In view of the rather delicate dependence of CQNS on these correlations (as discussed above), a realistic treatment of the kinematics in the hadron formation process is mandatory, including energymomentum conservation, noncollinear kinematics, and a well defined equilibrium limit. The recombination approach developed in Ref. (12) satisfies these requirements. In addition to this, the QGP evolution and subsequent hadronization are linked via the ansatz that resonances play an essential role in hot QCD matter around . This scenario is consistent with effective potential models where a nonperturbative description of the strongly coupled QGP (sQGP) is realized via bound (23) and/or resonance (24); (25) states of deconfined partons. Recent lattice QCD computations support the picture of various (light and strange) hadronic states surviving up to temperatures of 1.52 (26); (27).
Our article is organized as follows. In Sec. II we review our earlier developed model (12) for resonance hadronization based on the Boltzmann equation. In Sec. III we elaborate the computation of the phasespace distributions of quarks obtained from Langevin simulations of an expanding QGP fireball at RHIC. In Sec. IV we discuss the numerical results for the coefficients of and mesons within our model, and discuss their properties in terms of CQNS in both transverse momentum and transverse kinetic energy, . Sec. V contains our conclusions.
Ii Recombination from the Boltzmann Equation
Following Ref. (12), our description of hadronization at the critical temperature, , is based on the Boltzmann equation using resonance quarkantiquark cross sections to compute meson spectra in terms of underlying anti/quark phasespace distributions, (baryons could be treated in a similar way, e.g., in a twostep process using subsequent quarkquark and quarkdiquark interactions; in the present paper, we will focus on mesons). The meson phasespace distribution, , is determined by the equation
(1) 
where and denote threemomentum and position of the meson, , and (, =: meson mass and energy). The total meson width, , is assumed to be saturated by the coupling to quarkantiquark states, and taken to be constant, with the factor accounting for Lorentz time dilation, see also Ref. (11). Integrating over the fireball volume leads to the momentumdistribution function of the meson, , and the pertinent transport equation
(2) 
The drift term vanishes upon integration over since it is a total divergence: . The relation of the gain term, , to the underlying microscopic interaction is given by
(3) 
with the cross section for the process at centerofmass (CM) energy squared, , where are the fourmomenta of quark and antiquark. The quark phasespace distribution functions are normalized as . Throughout this paper, quarks will be assumed to be zerowidth quasiparticles with an effective mass (which contains both thermal and bare contributions). The classical nature of the Boltzmann equation warrants the use of classical distribution functions for all the particles, and we assume zero chemical potentials for all quark species. The cross section is approximated by a relativistic BreitWigner form,
(4) 
where is a statistical weight given in terms of the spin (color) degeneracy, (), of the meson (anti/quark), and denotes the quark threemomentum in the CM frame. With being the only channel, it follows that . Detailed balance requires the same in the loss term on the righthand side of Eq. (1), thus ensuring the correct equilibrium limit with the pertinent relaxation time. This formulation conserves fourmomentum and applies to all resonances with masses above the threshold, i.e., for a positive value,
(5) 
If the channel proceeds too far offshell, i.e., and (e.g., for pions), other processes need to be considered, e.g., (which, in principle, is possible in the present framework by implementing the respective cross sections). We note that the majority of the observed pions are believed to emanate from resonance decays (, , etc.); in addition, hydrodynamic calculations suggest that the elliptic flow in heavyion collisions at RHIC does not change much after hadronization (28). We also note that in the absence of a confining interaction individual quarks and antiquarks remain a part of the heat bath.
The equilibrium limit is readily recovered within our approach by imposing the stationarity condition,
(6) 
Then (2) is immediately solved by
(7) 
which represents the large time limit of the Boltzmann equation and is the expression that comes closest to the conventional QCM approximation. For hadronization times less or comparable to the relaxation time, , the equilibrium limit will not be reached and an explicitly timedependent solution is in order. We have verified numerically that Eq. (7) accurately recovers the standard thermal Boltzmann distribution for a meson at temperature , if the constraint of a positive value is satisfied (for negative the channel is inoperative), as illustrated in Fig. 1 for the case of the meson: when using thermal quark input distributions with radial flow (blastwave model), the computed equilibrium expression (7) is in excellent agreement with the same blastwave expression (identical temperature and flow profile) directly applied at the meson level. This reiterates the close connection between equilibration and energy conservation in the approach of Ref. (12), providing a significant improvement of previous QCMs.
Iii Partonic Spectrum
While the blastwave model provides a convenient description of the thermal component of empirical hadron (and/or parton) spectra, the transition to the kinetic, and eventually hardscattering, regime is more involved, especially with regard to phasespace correlations within QCMs, as discussed in the Introduction. In an attempt to generate realistic quark input distributions for meson formation processes in the vicinity of , we here adopt a FokkerPlanck approach for test particles evolving in a thermally expanding QGP background. The latter is parameterized with guidance from hydrodynamic models for central and semicentral AuAu collisions at RHIC, implementing empirical values (i.e., adjusted to experiment) for bulk matter properties such as total entropy, radial and elliptic flow. The elliptic flow is parameterized in terms of a flow profile, whose direction at each transverse position is perpendicular to confocal elliptic isobars. Its magnitude is chosen to increase linearly in distance from the center with an average boundary value of at the end of the mixed phase of the fireball evolution, i.e., at the hadronization time. The acceleration is adjusted so that the fireball at hadronization is approximately circular, with bulk . This model has been employed before in the context of heavyquark (HQ), i.e., charm and bottom spectra at RHIC (21), and good agreement with hydrodynamic simulations has been found (22) for the same HQ diffusion coefficient. HQ observables, which at present are semileptonic singleelectron decay spectra (29); (30), exhibit an unexpectedly large suppression and elliptic flow which cannot be understood within perturbative QCD (pQCD) including both radiative and elastic scattering. The key microscopic ingredient in Ref. (21) are resonant heavylight quark interactions mediated via effective (broad) and mesons in the QGP (31), inspired by the findings of thermal lattice QCD. In connection with a “conventional” coalescence afterburner (32) at , the predictions for singleelectron suppression and turned out to be in fair agreement with data (29); (30). In more recent work (25), the effective interactions have been replaced by inmedium matrices based on finitetemperature HQ potentials extracted from lattice QCD; these calculations not only confirmed the interaction strength generated by heavylight resonances in an essentially parameterfree way, but also identified prehadronic meson and diquark channels as the most relevant ones. This, in turn, provides a direct link between two main discoveries at RHIC, namely the strongly interacting nature of the sQGP and quark coalescence from a collective partonic source.
In the present paper, we build upon the above findings by extending the
FokkerPlanck approach to strange () quarks. While its applicability
criterion, (: transverse mass,
: momentum transfer in a typical scattering), seems to be only
marginally satisfied for momenta (at least in the early
phases of the QGP evolution), we note that most of the fireball
evolution occurs for temperatures close to . At higher (or
) one enters the kinetic regime where the FokkerPlanck treatment
becomes reliable again. With these limits properly satisfied, one may
hope that the Langevin simulations also accomplish a reasonable
description of the low regime () for
strange quarks, including realistic spacemomentum correlations, which
is one of the main objectives in our work. It remains to specify the
interaction strength of the strange and charmquark species. Here we
take further guidance from phenomenology by requiring that the final
quark exhibits the characteristic saturation (or maximum)
value of . Our baseline interaction for the stochastic
Langevin force is elastic pQCD scattering with a rather large value of
(which can be thought of as containing radiative and/or
parts of nonperturbative contributions). However, additional
nonperturbative interactions are necessary to achieve a sufficiently
large . As in Refs. (31); (21) we
associate these with mesonic resonance states with an interaction
strength controlled by the resonance width (with a larger width implying
stronger coupling). For quarks () the “heavylight”
resonances require a width of , and
for quarks (), which
is compatible with
Refs. (31); (21)
The quark phasespace distributions resulting from the Langevin approach embody strong correlations between spatial and momentum variables. E.g., quarks at high tend to be located in the outer layers of the fireball with a preferential alignment of the momentum and position vector directions (i.e., the quark momentum tends to point “outward”). Likewise, the collective (radial and elliptic) flow, which implies a welldefined (hydrolike) correlation between the position of the fluid cell and its radial motion, imprints this correlation on the (partially) thermalized components of the Langevingenerated quark spectra. We recall again that the often used factorized implementation of the coefficient in coalescence models completely ignores these rather elementary dependencies.
The proper implementation of the differential phasespace information, carried by the quark distribution functions (which is essential for a realistic discussion of hadronic scaling properties), into the hadronization formalism requires a few technical remarks. Since thermalization in the longitudinal direction is somewhat controversial, we assume the quark distributions to be homogeneous in the spatial coordinate and flat in rapidity. This leaves four independent transverse variables for each particle, which we choose in azimuthal form, , corresponding to the distribution . This 4D phasespace is then divided into finite bins (with a maximum value of ). For each simulated test quark, its final location and momentum is sorted into this grid. To warrant a (statistically) sufficiently smooth behavior of the computed meson observables, a sample of test particles is needed. Finally, an interpolation algorithm has been devised for converting the discretized distribution back into a continuous function, to be plugged into the hadronization formula. The algorithm recovers the periodicity properties of the two angular variables, and converges to an arbitrary sampled function in the limit of a large number of grid points. A suitable grid dimension corresponding to the above variables amounts to, e.g., , where the large number of points in is dictated by the large sensitivity in the determination of the elliptic flow coefficient, .
We furthermore have to specify how to treat partons that escape the fireball prematurely, i.e., before the end of the QGP/mixed phase is reached. Clearly, these partons preferentially carry a high , and, upon exiting the fireball, could undergo a hadronization mechanism different from coalescence, such as fragmentation. Since our fireball is isotropic, the transition from the QGP to the vacuum is a sharp one and it would be unrealistic to coalesce the exiting quark with a thermal distribution at a temperature above , or at a time much later than the exit time (when the fireball has cooled down to ). We therefore decide to include only partons in our hadronization framework which remain inside the fireball throughout the entire QGP evolution. This leads to an underestimation of the highmomentum part of the hadronic spectra. A comprehensive calculation for quantitative comparison to experiment also at high would require the treatment of the exiting partons (hadronized with either fragmentation or coalescence). Similarly, the hadronic we compute only reflects the partons within the QGP fireball at the end of its lifetime (which, however, do include nonthermal components from the Langevin simulation, in addition to the thermalized part of the spectrum).
The generated spectrum, which originally represents a probability distribution, requires a suitable normalization. Since the empirical light and strange hadron spectra are consistent with chemical equilibrium close to the expected phase boundary, we assume this to apply at the quark level as well, at the critical temperature of our fireball evolution. The fireball volume has been adjusted to match the total entropy of the fireball to the experimental hadronic final state multiplicities at at given collision centrality, e.g., for semicentral AuAu collisions (note that one fireball covers approximately units in rapidity). For charm quarks we augment the chemical equilibrium number by a fugacity factor, (33), to match their number to the expected hard production in primordial nucleonnucleon collisions (binary collision scaling; in Ref. (12) at leads to the same number of pairs in the fireball). We note, however, that the overall normalization has little impact on our main considerations of spacemomentum correlations and systematics.
Let us finally specify the assumptions on the rapidity () distributions. As mentioned above, for both charm and strange quarks we employ a step function as
(8) 
where the parameter depends on the quark mass and has been adjusted to recover approximately the fullwidthhalfmaximum of a thermal spectrum (amounting to and in connection with the quark masses quoted above).
Fig. 2 summarizes the results of the Langevin simulations for the probability distributions (integrated over spatial coordinates) for and quarks, compared to (normalized) blastwave spectra for , using the same flow field as in the Langevin simulation. The average surfaceexpansion velocity is . At low the spectra approach the equilibrium limit as to be expected since the background medium for the Langevin simulation is assumed to be fully equilibrated at all , as in a hydrodynamic calculation, and the stochastic process has been realized as to guarantee the correct equilibrium limit (including the adjustment of the (longitudinal) diffusion coefficient according to Einstein’s fluctuationdissipation relation (21)). As discussed above, the interaction strength has been chosen to recover the empirically observed maximum elliptic flow of . Note, however, that the assumption of a fully thermalized background medium implies rather large values at high , while the phasespace density of thermal partons is rather small. A full treatment of this problem would require to solve a selfconsistency problem, where the parton spectra in the background fireball evolution also exhibit a saturation of the elliptic flow, . This is beyond the scope of the present paper.
Iv Meson Spectra and Systematics
We now combine the ingredients of our resonance recombination model (RRM) by implementing the quark spectra computed from Langevin simulations in the previous Section (III) with the Boltzmannbased hadronization formalism of Sec. II, evaluated in the stationary (equilibrium) limit according to Eq. (7). The key issue we address is how the properties of the input (nonobservable) quark spectra reflect themselves in the (observable) meson spectra, in particular whether the spacemomentum correlations generated in the Langevin simulations can be consistent with the empirically observed CQNS, which, in turn, opens a window on the quark spectra at hadronization.
It remains to specify the masses and widths of hadrons in the recombination process. In line with our restriction to mesons located above the quarkantiquark threshold (due to the limitation to processes) we consider ,  coalescence with resonance masses corresponding to the vacuum values for () and () mesons; in connection with the quark masses as given above this implies similar values of . The (total) meson widths are chosen of comparable magnitude, i.e., and . As elaborated in Ref. (12), the numerical results, especially for the meson , are rather insensitive to variations in the meson width as long as is positive and substantially smaller (not smaller) than the resonance mass (width).
In Fig. 3 we display our RRM results for spectra of and mesons, including available RHIC data. Overall, the spectra largely agree with those computed in our previous work (Fig. 1 in Ref. (12)), where the input spectra were solely based on a blastwave parameterization. This is not surprising, since the quark spectra employed in the present work show a rather large degree of thermalization, up to momenta of  for strange and charm quarks, cf. Fig. 2. Consequently, the meson spectra shown here are somewhat harder than in Ref. (12) beyond due to the presence of the kinetic (hard) components in the quark spectra resulting from the Langevin evolution. The computed spectra for the are quite reminiscent to earlier blastwave based results (32); (34); (35). Note, however, that in Ref. (35) the recombination (blast wave) contribution only amounts to about of the total yield, significantly less than in the present paper (this is sensitive to the total opencharm cross section, which is not very well determined yet).
The RRM results for the meson are summarized in Fig. 4 (solid lines) and compared to the underlying quark scaled to meson variables in the conventional (empirical) way as . We find that for both the and the agreement is rather impressive, within a few percent relative deviation. The wiggles at the quark level are, to a large extent, driven by the finite grid sampling due to a step width of (the statistical error inherent to the Langevin simulation is smaller than that, using test particles). The convolution of two quark distributions results in much smoother curves at the meson level (due to the fitting procedure of the quark input). We are thus able to approximately recover CQNS in a microscopic calculation with the full information on spacemomentum correlations, characteristic for hydrodynamic expansion at low and a kinetic regime at higher .
Another potential source for scaling violations are flavor (or mass) dependencies at the quark level (rather than in the coalescence process). In Fig. 5 we compare the elliptic flow of different quarks (left panel) and mesons (right panel) with each other. The left panel confirms that the quark deviates significantly (up to at low ) from the strangequark . At high , the strangequark is a bit high, which, in principle, could be readjusted by a somewhat reduced strength of the resonance interactions in the Langevin simulations. At the meson level, the differences are similar. We have verified that comparable deviations persist when reducing the strange quark .
Finally, we address the question of CQNS with respect to transverse kinetic energy, , rather than transverse momentum, . Such a scaling has recently been highlighted by the PHENIX (19) and STAR collaborations (20) and seems to be very well satisfied by all available RHIC data, in centrality, collision energy and collision systems, after a geometric correction for the nuclear overlap. scaling has also been found to result from certain classes of hydrodynamic solutions (40), and therefore been argued to reflect a collectively expanding thermalized system of partons (41). From the point of view of quark coalescence, the problem of reconciling quark distribution functions with spacemomentum correlations (as implied by hydrodynamic expansion) with CQNS persists. In Fig. 6 we display the RRM results for and for the two different flavors as a function of . We find that the quark input distributions from the QGP Langevin simulations indeed share a rather universal behavior up to , encompassing both the quasiequilibrium regime at low energies and the kinetic regime at intermediate energies characterized by a leveling off at , cf. left panel of Fig. 6. We recall that the only adjusted input to this result is the common maximum value of the individual quark elliptic flow at about  (as suggested by the empirical CQNS deduced from experiment), controlled by the nonperturbative interaction strength in the stochastic Langevin force (again, a fine tuning for the quark would improve the agreement at higher ). The approximate universality at the quark level is nicely preserved at the meson level as a result of our Boltzmannbased recombination formalism, see right panel of Fig. 6. This is a quite remarkable result in view of the underlying spacemomentum correlations in our approach, which has not been achieved before in this form.
V Conclusions
In the present paper we have extended our previously formulated quark coalescence formalism, utilizing resonance interactions within a Boltzmann equation, by implementing microscopic quark phasespace distributions generated via Langevin simulations of an expanding thermal QGP fireball for AuAu collisions at RHIC. In this way we could combine the merits of our recombination approach (energy conservation and a proper equilibrium limit) with those of realistic quark distributions, which in particular encode the transition from a thermal regime at low to a kinetic one at intermediate . The latter feature is especially important as it produces the levelingoff of the elliptic flow, a key feature of observed hadron spectra and a crucial prerequisite to test any kind of quark scaling behavior. The (constituent) quarkmass and meson parameters were fixed at rather standard values, and we have constrained ourselves to mesons which are reliably calculable in our recombination setup (i.e., for positive values). The only real adjustment concerned the interaction strength in the Langevin process, as to reproduce the empirical maximum value for the quark (with the QGP fireball parameters tuned to empirical values of radial and elliptic flow, as in earlier applications to, e.g., heavyquark observables). These interactions have been modeled via (meson) resonances in the QGP, which we identified with the states formed in the coalescence process at , leading to the notion of a “Resonance Recombination Model” (RRM). Since the FokkerPlanck approach as an expansion of the Boltzmann equation is strictly valid for sufficiently massive and/or highmomentum particles, we restricted ourselves to “heavy” flavors, i.e., charm and strange quarks. At low , the latter are at the borderline of applicability of a FokkerPlanck framework.
Our main finding is that within this rather generic setup, largely based on first principles augmented by a concrete realization of the strongly interacting QGP, the constituent quark scaling of the meson elliptic flow emerges rather naturally including spacemomentum correlations characteristic for a collectively expanding source. The scaling holds for individual mesons, but appears to be rather universal in quark and meson flavor (mass), especially when applied in (transverse) kinetic energy rather than momentum, which is in line with recent experimental findings. By overcoming some of the limitations of previous (more schematic) coalescence models, and by achieving the first robust implementation of realistic (microscopically computed) phasespace distribution functions of quarks, our formalism could provide a useful tool to better understand systematics of RHIC data, most notably the interplay of a thermal and kinetic regime in connection with phasespace properties of the partonic fireball as viewed through the hadronization process.
Clearly, a formidable list of open issues persists, including the extension to light quarks, the role of gluons and of deeply bound hadronic states (possibly requiring additional formation processes), more realistic spectral functions of mesons and quarks and their interactions (both around and above), a selfconsistent treatment of thermal and kinetic components (possibly requiring full parton transport), a systematic classification of viable parton phasespace distributions, hadronic reinteractions, etc. Progress has already been made on a number of these aspects, but a comprehensive approach remains a challenging task.
Acknowledgements.
We thank T. Hahn for insightful clarifications about the CUBA multidimensional integration package (42) which has been used in this work, and R. J. Fries for valuable discussions. This work was supported in part by a U.S. National Science Foundation CAREER award under grant no. PHY0449489 and by the A.v.Humboldt Foundation (through a Bessel Research Award).Footnotes
 We recall (12) that the only requirement on the quark and meson masses is that the latter are above the twoquark threshold; within this restriction variations in the mass and (positive) values have little effect on the recombination process. In fact, as we will see below, CQNS scaling emerges approximately independent of quark mass.
References

K. P. Das and R. C. Hwa, Phys. Lett. B 68, 459 (1977), erratum ibid. 73, 504 (1978).
 E. Braaten, Y. Jia, and T. Mehen, Phys. Rev. Lett. 89, 122002 (2002).
 R. Rapp and E. V. Shuryak, Phys. Rev. D 67, 074036 (2003).
 T. S. Biro, P. Levai, and J. Zimanyi, Phys. Lett. B 347, 6 (1995).
 R. C. Hwa and C. B. Yang, Phys. Rev. C 67, 034902 (2003).
 V. Greco, C. M. Ko, and P. Levai, Phys. Rev. C 68, 034904 (2003).
 R. J. Fries, B. Müller, C. Nonaka, and S. A. Bass, Phys. Rev. C 68, 044902 (2003).
 D. Molnar and S. A. Voloshin, Phys. Rev. Lett. 91, 092301 (2003).
 Z.W. Lin, C. M. Ko, B.A. Li, B. Zhang, and S. Pal, Phys. Rev. C 72, 064901 (2005).
 J. Zimanyi, P. Levai, and T. S. Biro, J. Phys. G 31, 711 (2005).
 H. Miao, C. Gao, and P. Zhuang, Phys. Rev. C 76, 014907 (2007).
 L. Ravagli and R. Rapp, Phys. Lett. B 655, 126 (2007).
 A. Ayala, M. Martinez, G. Paic, and G. T. Sanchez, Phys. Rev. C 77, 044901 (2008), eprint 0710.3629.
 W. Cassing and E. L. Bratkovskaya, Phys. Rev. C 78, 034919 (2008).
 R. J. Fries, V. Greco, and P. Sorensen, Ann. Rev. Nucl. Part. Sci. 58, 177 (2008).
 S. Pratt and S. Pal, Nucl. Phys. A 749, 268 (2005).
 D. Molnar (2004), eprint arXiv:nuclth/0408044.
 V. Greco and C. M. Ko (2005), eprint arXiv:nuclth/0505061.
 A. Adare et al. (PHENIX Collaboration), Phys. Rev. Lett. 98, 162301 (2007).
 B. I. Abelev et al. (STAR Collaboration), Phys. Rev. C 75, 054906 (2007).
 H. van Hees, V. Greco, and R. Rapp, Phys. Rev. C 73, 034913 (2006).
 G. D. Moore and D. Teaney, Phys. Rev. C 71, 064904 (2005).
 E. V. Shuryak and I. Zahed, Phys. Rev. D 70, 054507 (2004).

M. Mannarelli and R. Rapp, Phys. Rev. C 72, 064905 (2005).
 H. van Hees, M. Mannarelli, V. Greco, and R. Rapp, Phys. Rev. Lett. 100, 192301 (2008).
 F. Karsch and E. Laermann (2003), published in R. C. Hwa, X. N. Wang (editors), Quarkgluon plasma vol. 3 (World Scientific, 2004) p1, eprint arXiv:heplat/0305025.
 M. Asakawa and T. Hatsuda, Nucl. Phys. A 721, 869 (2003).
 P. F. Kolb and U. W. Heinz (2003), published in R. C. Hwa, X. N. Wang (editors), Quarkgluon plasma vol. 3 (World Scientific, 2004) p634, eprint arXiv:nuclth/0305084.
 A. Adare et al. (PHENIX Collaboration), Phys. Rev. Lett. 98, 172301 (2007).

B. I. Abelev et al. (STAR Collaboration), Phys. Rev. Lett. 98, 192301
(2007).
 H. van Hees and R. Rapp, Phys. Rev. C 71, 034907 (2005).
 V. Greco, C. M. Ko, and R. Rapp, Phys. Lett. B 595, 202 (2004).
 L. Grandchamp, R. Rapp, and G. E. Brown, Phys. Rev. Lett. 92, 212301 (2004).
 A. Andronic, P. BraunMunzinger, K. Redlich, and J. Stachel, Nucl. Phys. A 789, 334 (2007).
 X. Zhao and R. Rapp, Phys. Lett. B 664, 253 (2008).
 A. Adare et al. (PHENIX Collaboration), Phys. Rev. Lett. 98, 232301 (2007).

S. S. Adler et al. (PHENIX Collaboration), Phys. Rev. C 72, 014903
(2005).
 J. Adams et al. (STAR Collaboration), Phys. Lett. B 612, 181 (2005).
 B. I. Abelev et al. (STAR Collaboration), Phys. Rev. Lett. 99, 112301 (2007).
 M. Csanad et al., Eur. Phys. J. A 38, 363 (2008).
 R. A. Lacey and A. Taranenko, PoS C FRNC2006, 021 (2006).
 T. Hahn, Comput. Phys. Commun. 168, 78 (2005).