Quarkonium production and suppression: Theory

Heavy quarkonium theory has seen significant progress over the past years. Not only has the community steadily improved the understanding of fully equilibrated quarkonium using first principles lattice QCD studies, but it has taken major steps towards developing a genuine dynamical understanding of quark-antiquark pairs immersed in a hot medium by embracing the paradigm of open quantum systems. We review recent developments, placing an emphasis on quarkonium real-time dynamics.


Medium effects
. Overview of the different physics mechanism at play for in-medium heavy quarkonium that perturbation theory is applicable (nb. non-perturbative nuclear parton distribution functions however play an important role here too). In addition QQ production factorizes from the quarkonium formation process, which itself is genuinely non-perturbative. Theory is making progress in quantitatively elucidating quarkonium production in p+p and p+A collisions by improving on phenomenological models, such as the color evaporation model [10] and by recently marrying the effective field theory NRQCD and the color glass condensate framework [11]. A similarly robust understanding of early production and formation in A+A has however not yet been established. If a quarkonium state forms early on, it finds itself soon immersed in the quark-gluon-plasma formed by the light remnants of the collision. The incessant kicks of this hot environment will interfere with the binding of the quark-antiquark pair, which may either survive or dissociate along the way. Once the bulk matter reaches the crossover transition to the hadronic phase, colored partons have to combine into color neutral hadrons and we can entertain two possible scenarios: either a surviving in-medium state transitions into a vacuum state or (as we now know is relevant for charmonium at LHC) if a large enough number of quark-antiquark pairs was produced early on, there exists a finite probability for them to recombine into quarkonium even if they became fully dissociated in the QGP. Theory is actively exploring the consequences of a scenario in which all quarkonium is produced at the phase transition based on the statistical model of hadronization [12] and at this conferences new ideas on the coalescence of single heavy quarks into open heavy flavor during hadronization have been presented [13,14].
In the following let me concentrate on the recent progress made in understanding heavy quarkonium in the intermediate stages of the collision, where one considers a fully formed state ineracting strongly with its hot environment. To get an intuition for the physics involved, I invite the reader to think of heavyquarkonium in terms of an energy level diagram, often used in atomic physics and as shown on the left of Fig.1. Already Matsui and Satz [15] knew that in the presence of a medium the Hamiltonian of the quarkonium system changes, i.e. some former bound stationary states become unbound and the energy of the remaining ones may change (blue to red levels). What has been understood however over the past 12 years is that the medium also acts as a radiation field (orange arrow) which may introduce absorption and emission processes. Imagine we start with an color-singlet eigenstate of the in-medium Hamiltonian (orange dot, center box). Instead of remaining set, the medium will induce transitions to higher or lower lying states, as well as between the singlet and octet sector. In thermal equilibrium the transitions among states are balanced out, such that the probability to find the particle at a certain energy remains time-independent. Theory has non-perturbative access to this physics via the study of in-medium equilibrium spectral functions ρ(ω) from e.g. lattice QCD. At T = 0 they contain delta-like peaks where bound states exist (plus a continuum) while at T > 0 the peaks broaden and shift. The latter indicates the change in the Hamiltonian, the former encodes the probability to transition away from the state. In a heavy-ion collision the situation is different, as one starts out with an unequilibrated state (blue dot, top left figure). Now both the Hamiltonian induces mixing between the states, as well as the radiation field and these transitions are not balanced. Only at late times will the system equilibrate to time independent occupations. Here we are interested in the real-time evolution of the occupation numbers, or what is conventionally termed survival probabilities of different states. A lot of progress has been made to capture this physics, but a fully non-perturbative approach is still outstanding.
Extracting equilibrium spectral functions from lattice QCD is challenging, as it amounts to a highly illposed inverse problem. A lattice QCD simulation acts similarly to an imperfect detector producing a spectrum that is folded with a resolution function. Actually lattice QCD encodes ρ(ω) in so-called Euclideantime correlation functions D(τ), which are essentially their Laplace transform. Two strategies exist on the market to extract the spectral function: on the one hand we may use Bayesian inference, a data analysis technique which exploits prior QCD domain knowledge (e.g. positivity of ρ(ω)) to regularize the inverse Laplace transform. The well known Maximum Entropy Method (MEM) [16] or the Bayesian Reconstruction (BR) method [17] are two different incarnations of this strategy. While they do not presuppose a particular form of the spectral function they require very high precision input data to succeed. On the other hand one may model the spectral function based on prior information obtained from EFT's and perturbation theory and then fit the model parameters to the lattice data. This strategy gives access to e.g. excited state and transport physics but introduces model dependencies.
The Bayesian inference community has made lot's of progress understanding the systematic effects of the reconstruction process over the past years and is converging on consistent results among different groups e.g. for the stability estimates of ground state quarkonium [18,19]. While recent studies based on non-relativistic lattice EFT's have provided consistent results indicating that the mass shifts of quarkonium are negative, it remains an active research question. Recent results on T > 0 quarkonium dispersion relations were presented in [20]. All Bayesian studies so far agree that in-medium modification occurs in a sequential fashion, related to the vacuum binding energy of the quarkonium state. It needs to be kept in mind that this of course does not directly translate into sequential suppression in a heavy-ion collision. Limited by simulation data precision, so far no results on excited states and only limited information of spectral widths has been obtained.
Two quarkonium spectral modelling studies in lattice QCD have been presented at this conference. In [21] a non-relativistic formulation of heavy quarks is deployed on realistic simulations (m π ≈ 160MeV) of the thermal QCD medium background. Adopting techniques from T = 0 simulations, i.e. the use of extended sources for a better signal to noise ratio, to T > 0 and utilizing a simple parametrization of the spectral function, the authors fit the lattice data to obtain first results on the in-medium modification of bottomonium excited states and the width of the ground state. Again hierarchical ordering with vacuum binding energy is observed. The second study [8] considers fully continuum extrapolated simulation data based on a relativistic formulation of heavy quarks, which however are immersed in a purely gluonic medium only. The spectral functions here are modelled using combined prior information from perturbation theory and the EFT pNRQCD. Fits of the lattice data reveal a negative mass shift and allow to give a robust lower limit on the transport coefficient related to heavy quark diffusion D > 2/2πT .
Besides obtaining spectral information from lattice QCD, the definition and extraction of an intuitive in-medium heavy quark potential has been one focus within heavy quarkonium theory. We now understand how to define it systematically from QCD by approximating the physics of heavy quarks via a set of EFT's called NRQCD and pNRQCD (see e.g [22]). In the latter, quarkonium is represented as color singlet and octet wavefunctions and the Lagrangian of these degrees of freedom directly reflects the physical processes discussed in the center panel of Fig.1. The in-medium energy levels are given by a Schrödinger-like part in the Lagrangian, where the singlet and octet evolve according to their respective in-medium potentials. In addition radiation field effects of the gluons are included by dipole-like terms, which may induce transitions between the two sectors.
The parameters of this EFT, the interactions potentials and the strength of the dipole terms are set via a matching procedure, where we compute a correlation function in the EFT and a correlation function with the same physics content in QCD and set them equal at the energy scale of interest. In the limit of static quarks the so-called real-time Wilson loop is the QCD counterpart to the pNRQCD correlation function of two color-singlet wavefunctions (see e.g. [4,23,24]). Thus their real-time evolution can be identified and used to unambiguously define what we mean by the in-medium heavy quark potential in the language of EFTs. Note that this quantity in general can be complex. The real-time definition cannot directly be evaluated in lattice QCD simulations, which are carried out in an artificial Euclidean time but strategies have been developed [5,6,25], based on spectral function reconstruction, to extract its values non-perturbatively. The most recent results on the real-and imaginary part of the potential, based on Bayesian spectral reconstructions and a Pade approximation technique were given in [7]. This real-time potential has also been computed using methods of holography, presented as poster at this conference (see e.g. [26]).
Other definitions of potentials exist in the literature, which are not based on the above mentioned EFT's. E.g. the T-matrix approach [27] sets out to describe the physics of both the light and heavy quark degrees of freedom with a model Hamiltonian, which includes a kinetic term, as well as a four-fermi interaction, whose strength is given by a real-valued in-medium potential. The functional form of this potential is fitted, such that the model reproduces thermodynamic properties of QCD computed on the lattice, e.g. the QCD free energies. Remarkable agreement with the non-perturbative data is found. Since however the T-matrix potential is not based on a separation of scales and attempts to summarize the the physics of light and heavy quarks it is conceptually different from the EFT potential and thus it is neither surprising nor of concern that the best fit values for the T-matrix potential differ from the lattice EFT potential.
The study of the heavy quark potential in lattice QCD is one way how to build a first-principles bridge from QCD at finite temperature to the real-time dynamics of heavy quarkonium. First-principles insight is indeed needed, as the currently available experimental data (mostly R AA of quarkonium ground state particles) can be reproduced by a variety of phenomenological models, which incorporate widely different physics mechanisms. We can distinguish two large classes of studies: one uses a rate equation where a temperature dependent decay rate drives the system towards an equilibrium distribution (see e.g. [28][29][30]). The other, deployed mainly for bottomonium, solves a deterministic Schrödinger equation in the presence of a complex in-medium potential (see e.g. [31,32]). If we wish to quantitatively infer the physics underlying the experimental results and distinguish between different models, we need to establish their actual range of validity and interrogate QCD about appropriate parameters to use in these models. In addition, theory has the goal to derive more and more accurate real-time descriptions of quarkonium directly from QCD.
An important step towards establishing a general real-time approach to heavy quarkonium has been achieved in the community by embracing the open quantum systems framework. Originally developed in the context of condensed matter theory, it focusses on scenarios in which a small subsystem is immersed in a large (thermal) environment. The overall system, consisting of environment and subsystem, is closed. It possesses a hermitean Hamiltonian, which itself is composed of a part describing the interactions among the quark-antiquark pair H QQ , another part for the medium degrees of freedom H med and the interaction part between the two (the latter is often written as H int = m Σ m ⊗ Ξ m , the former acting in the quarkonium subspace, the latter in the medium subspace). The overall density matrix ρ(t) evolves according to the von-Neumann equation What we are interested in however is only the dynamics of the small QQ subsystem. Formally tracing out the medium degrees of freedom leads to the so-called reduced density matrix ρ QQ = Tr med [ρ]. The goal of theory is to derive from QCD the equation of motion for ρ QQ , called the master equation. A generic feature is that tracing out the medium leads to dissipative dynamics of the subsystem, even if the evolution of the overall system is fully unitary. The form of these master equations depends on the particular separation between three characteristic time-scales: The environment relaxation scale τ E from Σ m (t)Σ m (0) ∼ e −t/τ E , the intrinsic QQ scale τ S from τ S ∼ 1/|ω − ω | and the QQ relaxation scale τ rel , which for a heavy particle can be defined from the isotropization timescale of its momentum p(t) ∼ e −t/τ rel . In case that the environment relaxes much faster than the quarkonium system, the time evolution can be approximated as Markovian. This in turn allows it to be cast in the general form of a Lindblad equation which conserves hermiticity ρ † QQ = ρ QQ , the probability interpretation T r[ρ QQ ] = 1 and positivity n|ρQQ|n > 0, ∀n during the time evolution. The modified HamiltonianH QQ again encodes the energy level changes of Fig.1, while the Lindblad operators L k are intimately related to the radiation field effects of the medium.
As shown in Fig.2 two major branches are being explored in this line of study at the moment. On the one hand theorists use their experience with the EFT pNRQCD at high temperature to derive a Lindblad equation applicable for a weakly coupled medium. Assuming the time scale separation of the quantum optical limit τ S τ rel and performing a gradient expansion it has been shown how to reduce that Lindblad equation to  Historically one of the most common approach to bottomonium modeling is to resort to solving a deterministic Schrödinger equation which is governed by an in-medium potential. From our discussion in Section 4.2 we saw in Eq. (238) that such an ansatz amounts to a truncation of the dynamics, neglecting dissipative e ects, and in addition to an adiabatic approximation that averages over the fluctuations. From the numerical tests within the stochastic potential model, we learned that in such an adiabatic truncation the survival of quarkonium states may be underestimated if the e ects of decoherence are relevant. I.e. the RAA computed in that way may underestimate the correct physical value.
Combining the deterministic Schrödinger equation with an anisotropic hydrodynamic description of the medium background evolution, as well as a Glauber model based initial distribution of primordial states and the decays accounting for feeddown after hadronization, the Kent State University group has explored the suppression of bottomonium in heavy-ion collisions in detail. A full description of the framework is found in Ref. Let us start with Bottomonium at STAR, since there the separation of scales between the medium temperature created in the collision center and the bottom mass is most pronounced. In addition bottomonium is expected to act as a true non-equilibrium probe, without reaching any significant degree of kinetic thermalization. In such a scenario, one may speculate that decoherence is not yet very e ective and thus the adiabatic approximation may prove already satisfactory. In the left panel of Fig. 55  J _ RAA of the QGP phase or whether they are e ciently melted on the way. Two groups, one based at Texas A& M and the other at Tsinghua university have developed transport models based on the Boltzmann equation, respectively on a corresponding averaged rate equation. Both approaches have in common that they implement the possibility for dissociation of primordial charmonium particles, as well as the dynamical recombination of such states throughout the QGP evolution. Feed-down from excites state is included at the end of the medium evolution. Non-perturbative information on quarkonium dissociation is implemented in the former approach by computing their stability properties via the T-matrix approach or via melting temperatures in the latter. As we discussed in Section 3.2.1, the definition of the potential used in a Bethe-Salpeter based approach and its relation to the EFT based potential remain an active area of research and constitute one significant source of uncertatinty. On the other hand we have seen that the definition of a melting temperature in the presence of a thermal width is not uniquely defined and especially di cult when using direct lattice QCD results, which in turn contributes to the overall systmatic uncertatinty. The recombination probability in these models is constructed using arguments of detailed balance, which are expected to work well close to equilibrium but may require corrections at early times far from equilibrium. More details on the construction of loss and gain terms can be found in Ref.
[337]. The medium evolution is treated somewhat di erently among the two models. In the former the concept of entropy conservation in conjunction with the measured particle multiplicities is used to model a temperature profile within an isotropically expanding fireball. The latter model on the other hand computes a temperature profile directly based on 2 + 1 dimensional Bjorken expansion in the QGP phase. In each case a first order transition like regime is used to connect the QGP with a hadron resonance gas phase at lower temperatures.
Both models are able to describe the RAA in terms of centrality and transverse momentum at RHIC and LHC in an equally good fashion as shown for the example of its centrality dependence in the center panel of Fig. 52, see also Ref.
[300] (in non-central collisions the model of Ref. [336] seems to somewhat underestimate the actual RAA). For the pt dependence we select here for better readability in the right panel of Fig. 52 a single results from Ref. [337] compared to the measurements by the ALICE collaboration at˘sAA = 2.76TeV. As indicated by the two di erent black lines referring to the primordial and the regeneration component of the total RAA there are two di erent regimes present, which are smoothly connected. At large pt primordial charmonium appears to contribute a majority of the yield, while at small pt a significant fraction of produced J _ arises from recombination e ects. Similar behavior for the RAA dependence on centrality is observed: in non-central collisions primordial J _ dominates but for central collisions regeneration plays an almost equally important role. The e ects of regeneration at LHC are pronounced but also at RHIC the transport model computations indicate that dissociation of primordial charmonium alone cannot account for the observed patterns in RAA.
From considerations of the nuclear modification factor of the charmonium ground state J _ we have so far learned that dissociation and regeneration are leading to an intricate pattern of charmonium suppression requiring insight into all stages of the collision. While in the transport models a partial equilibration of the charmonium states occurs, that equilibration is assumed to be complete in the statistical model. The question thus remains is there a way how to distinguish between these two scenarios even though both reproduce the RAA well. This question becomes even more  Overview the current open-quantum-system based approaches to quarkonium dynamics. Note that they provide a systematic path from QCD to the most popular phenomenological models on the market through a series of well specified approximations, establishing robustly a range of validity.
the Boltzmann equation and in case of further time-scale separations even further to a rate equation. As presented at this conference [33] this result was achieved by expressing the operators Σ and Ξ of the interaction Hamiltonian by the quantities present in the weak-coupling pNRQCD Lagrangian and relating the distribution function of the Boltzmann equation to the color singlet matrix elements of the reduced density matrix [34]. The important point here is that this derivation provides consistent pNRQCD based expressions for the recombination and dissociation scattering kernels in the Boltzmann equation. Simplifying the description of the color octet distribution functions by treating dissociated quarks and antiquarks via a Langevin diffusion equation, this Boltzmann equation has been deployed to estimate quarkonium yields in [35]. In order to proceed towards a non-perturbative treatment, another group has focused on the subset of deep Coulombically bound quarkonium states immersed in a strongly coupled medium [9]. It can be shown that in this specific case the dynamics are governed by two transport coefficients, the heavy quark diffusion constant, as well as a potential correction term. As they can be expressed in terms of correlation functions of color electric fields, these transport coefficients are amenable to an evaluation in lattice QCD, a task whose completion is ongoing. As is the case with all pNRQCD based computations, we must keep in mind that it is formulated for a specific hierarchy of scales, which needs to be ascertained when deploying these methods in a phenomenological simulation concurrently evolving many different quarkonium states.
Along the second, lower branch in Fig.2 the community works on explicitly tracing out the medium degrees of freedom in the path integral language. This leads to the so-called Feynman-Vernon influence functional, which encodes all information about QQ medium interactions. In general it represents an extremely complicated expression and only through further approximations can it be tamed, giving way to a Markovian Lindblad equation. Interestingly it is possible to show that such master equations are equivalent to a non-linear stochastic Schrödinger equation. This has allowed for the first time to give a QCD derivation to such non-linear approaches (see e.g. [36]) previously introduced in the phenomenology community. Going to the recoilless limit (neglecting dissipation) the dynamics can be reduced to a linear Schrödinger equation, which recovers another phenomenological model, the stochastic potential of Ref. [37]. Finally going to the adiabatic limit one recovers the deterministic Schrödinger equation currently deployed in the study of Bottomonium states (see e.g. [31]) . Again we are able to give a systematic chain of approximations and thus clear ranges of validity for phenomenological models thanks to the open quantum systems approach.
Let us focus on one strategy to derive the master equation, which focuses on high temperatures. There, the medium is weakly coupled and thus a perturbative expansion may be deployed. Assuming the time-scale separation of the quantum-Brownian motion limit τ E τ S , which is often realized among quarkonium states, the following momentum and color dependent Lindblad operator has been derived in Ref. [38,39] L k,a = D(k) I emphasize this equation, since it embodies an interesting connection between quarkonium real-time dy-namics and the EFT based heavy quark potential discussed before. The two terms correspond to the physics of the quark and antiquark respectively (overall color neutrality of quarkonium requires opposite sign). Those terms with 1 can be attributed to fluctuations feeding energy into the subsystem, the momentum dependent terms on the other hand encode dissipative effects that return energy to the medium. Together they allow the quarkonium particle to eventually thermalize with its surroundings. Both phenomena are governed by a function christened D, which is intimately related to the imaginary part of the complex EFT potential On the other hand, the Debye screened real part of the EFT potential enters in the modified HamiltonianH QQ of the corresponding Lindblad equation (2). Hence this QCD based result tells us how the complex heavy quark potential governs quarkonium dynamics at high temperatures, which clearly goes beyond a simple Schrödinger equation. The question of how to extend this approach, currently limited to high temperatures, to a strongly coupled medium is the focus of ongoing research efforts.
There is another bit of insight hiding in plain sight: the function D(r) exhibits a well peaked structure around the origin r = 0, which can be identified as a new intrinsic scale in the system, the correlation length corr of the medium (see e.g. discussion in Ref. [40]). In the context of quarkonium dynamics it plays an important role, as it characterizes the effects of decoherence on the quarkonium wavefunction. If the corr is much larger than the extent of our quarkonium state, L k,a will induce essentially the same color rotations for the Q andQ leaving the bound state almost unaffected. Once corr becomes similar to the quarkonium size, the quark and antiquark receive different color rotations and the binding of the state is strongly impaired. Why did we say that corr is related to decoherence though? When we have a closer look at the effect of the function D, we find that it acts approximately as ρ QQ (x 1 , i.e. it induces the decay of the off-diagonal terms of the density matrix, which is exactly what decoherence does. I.e. it is understood that the real-time dynamics of quarkonium arise from an intricate interplay of screening (characterized by the Debye mass m D ) and decoherence ( characterized by corr ).
If decoherence is efficient, one expects that a semiclassical description will eventually become applicable. As investigated in detail in [41,42], if in addition the color degrees of freedom thermalize quickly, the non-Abelian nature of the theory does not play an important role and one may go over to a Langevin-like description of the heavy quarkonium states. One of the attractive features of such a description is the possibility to straight forwardly generalize it to many QQ pairs, allowing insight into recombination dynamics which are difficult to compute in a fully quantum setting.
One aspect of the quarkonium dynamics that two poster contributions to this conference ( [43,44]) elucidated is the effect of decoherence in the presence of explicitly implemented color degrees of freedom. To this end both groups restricted themselves to the recoilless limit (the former in 1d, the latter in 3d), which is expected to be applicable at early times in the evolution. In this approximation the dynamics is represented by an ensemble of wavefunctions (that carry a color index) which are governed by a linear Schrödinger equation with real-valued attractive color singlet and repulsive octet potential, as well as two noise terms in color space. By inspecting the color singlet and octet density matrices, it was found that while decoherence is leading to a diagonal structure for large values of the separation distances, for small distances, remnants of the initial bound state may persist for an appreciable amount of time. Taking the Wigner transform of such a density matrix, the structures around the origin translate into negative contributions of the Wigner distribution function. This in turn signals that the genuine quantum nature of quarkonium can remain relevant during the evolution in a hot medium and that if we wish to apply a semi-classical approximation its applicability needs to be established on a case-by-case basis.
Another recent study presented as poster contribution at this conference [45] took on the task to implement the fully dissipative dynamics of the Lindblad equation at high temperature discussed above. As is known in the condensed matter open quantum systems community, any Lindblad master equation can be implemented (the jargon term is unravelled) in terms of an ensemble of wavefunctions that evolve according to a non-linear stochastic Schrödinger equation. Using the unravelling technique of quantum state diffusion, the authors found that simulations, starting from different initial states (ground or first excited state), will at late enough times approach a universal steady state of non-vanishing occupation numbers (within the finite volume of a heavy-ion collision). This is encouraging, as both the deterministic Schrödinger equation with complex potential and the dissipationless stochastic potential model all lead to vanishing bound state contributions at late times. One may ask whether the steady state that is observed is a thermal state, which the authors answer in the affirmative by plotting the late-time occupancies vs. the mode energy. This revealed a clean exponential Boltzmann behavior, whose exponent excellently reproduced the environment temperature. The open quantum systems framework considered here is a versatile theoretical tool, as it provides a unified language to describe the evolution of a single heavy quark, a quarkonium pair and in principle also that of multiple QQ pairs, even if the computational cost is extremely high. As a result e.g. the function D(r), related to the imaginary part of the interquark potential in the medium also governs the dynamics of the single quark (see e.g. [46]). A similar connection between quark transport and evolution of quarkonium we discussed in the context of the heavy quark diffusion coefficient and the pNRQCD based master equation.
Heavy quarkonium theory has entered a new and exciting era, in which the dissipative real-time evolution of quark-antiquark pairs takes center stage. Combining the open quantum systems framework with the arsenal of high precision lattice QCD simulations available, the community is working on non-perturbative implementations of quarkonium dynamics that connect QCD to actual phenomenology and experimental data with unprecendented carity. Together with the activities of the theory community in better understanding the initial production process (not to forget the advances in the lattice QCD determination of parton distribution functions), a more and more comprehensive picture of quarkonium production and suppression will emerge over the coming years. The prospect of the measurement of excited states in the upcoming run3 at LHC promises access to much more discriminatory observables, which together with the newly found insight on the range of validity of phenomenological models, will help us to distill the relevant physical processes underlying the measured yields. Exciting news about first measurements in heavy-ion collisions of exotic quarkonium states, such as the X(3872) presented by CMS at this conference [47] furthermore opens up a new venue to learn about the microscopic composition of these states, which have been studied by theorists at T = 0 intensively.