A first look at Bottomonium melting via a stochastic potential

We investigate the phenomenon of Bottomonium melting in a thermal quark-gluon plasma using three-dimensional stochastic simulations based on the concept of open-quantum systems. In this non-relativistic framework, introduced in [Phys.Rev. D85 (2012) 105011], which makes close contact to the potentials derived in effective field theory, the $b\bar{b}$ system evolves unitarily under the incessant kicks by the constituents of the surrounding heat bath. In particular thermal fluctuations and the presence of a complex potential in the EFT are naturally related. An intricate interplay between state mixing and thermal excitations emerges as we show how non-thermal initial conditions of Bottomonium states evolve over time. We emphasize that the dynamics of these states gives us access to information beyond what is encoded in the thermal Bottomonium spectral functions. Assumptions underlying our approach and their limitations, as well as the refinements necessary to connect to experimental measurements under more realistic conditions are discussed.


Bottomonium as QGP thermometer
Elucidating the properties of strongly interacting matter shortly after the Big-Bang is the central reason behind scrutinizing the collisions of heavy-ions [2] at relativistic energies at the RHIC and the LHC collider facilities. Under the extremely large energy densities present in the very early universe as well as in the center of these man made collisions, the fundamental matter constituents, quarks and gluons, are not confined within hadrons but instead form a quark-gluon plasma (QGP) [3,4].
Measurements of e.g. the elliptic flow of light hadrons at RHIC [5,6] already helped to establish that a very low shear viscosity to entropy ratio η/s 1/4π [7] is a prime characteristic of this high temperature phase of quantum chromo dynamics (QCD), just above the deconfinement transition. The precise determination of such bulk properties of the fireball created in the collision center, ideally crosschecked by multiple independent methods, is steadily progressing but remains a challenging endeavor.
Another phenomenologically interesting property is the temperature reached in heavyion collisions. On the one hand electromagnetic probes, such as thermal photons emitted from the QGP have been utilized to obtain estimates both for T RHIC ( √ s N N = 200GeV) = 221 ± 19MeV and T LHC ( √ s N N = 2.76TeV) = 304 ± 51MeV [8,9], which lie well above the deconfinement crossover transition T ∼ 150MeV determined in numerical lattice QCD studies [10,11]. In this study on the other hand we focus solely on strongly interacting objects and how they can be used to infer medium properties. Describing their evolution hence requires us to deal with the real-time dynamics of quarks and gluons at finite temperature.
In general, to measure the temperature T of a physical system, we can resort to either a static or a dynamic strategy. The former, a thermometer in the everyday sense, requires us to bring into contact with the system of interest a test system, which after having equilibrated over time to a common temperature we decouple and inspect regarding its own thermodynamics properties, such as its density, free energy or entropy. The dynamical approach on the other hand observes the not necessarily thermal properties of a test system as it approaches equilibrium with its surroundings, just as the speed of melting a sugar cube in a cup of tea reveals the temperature of the beverage. While on macroscopic scales both approaches are often applicable, the physics at the subatomic realm, in particular within the strongly interacting QGP at the center of a heavy-ion collision, poses two important challenges to choosing an appropriate tactic.
First we need to find appropriate probes, whose behavior is sensitive to the surrounding medium, while their underlying structure remains largely unaffected. I.e. we require probes whose characteristic internal scales, such as their masses M , lie well above those of the medium partons. The more pronounced this separation is, the more straight forward their theoretical description becomes, as there exists a small parameter (dividing the scale of the medium by the scale of the probe), which can be used as an expansion parameter to construct effective approximate descriptions.
Secondly the system we wish to investigate only exists for a minute amount of time τ QGP ∼ 5 − 10fm, and has reached local thermal equilibrium only after τ →T QGP ∼ 1fm. As we require a clear separation of scales for our probes M T QGP , we assume here that their thermalization time, i.e. the time it takes for their number density to become Bose-Einstein (or Fermi-Dirac) distributed, to be comparatively long τ →T probe ∼ M T QGP τ →T QGP . Hence we will resort to a dynamic strategy for temperature measurements, as we anticipate that observing the thermal properties of our probes will be difficult within the lifetime of the QGP.
One popular candidate for an appropriate kind of probe are the bound states of heavy quarks and anti-quarks, so called heavy quarkonium. Here we will work with the Bottomonium family, in particular the vector mesons consisting of a b andb quark, with a constituent mass m b = 4.66 ± 0.03GeV (1S) [12].
The clear separation between m b T LHC , T RHIC underlies the intuitive notion that Bottomonium and its interactions with the surrounding medium can be captured in an effective non-relativistic description. More precisely, one wishes to summarize the in-medium physics of the two-body system in an effective potential V EFT (r) entering a non-relativistic Schroedinger equation. The formalism of effective field theory [13][14][15][16][17] adapted to finite temperature over the last decade, allows to put this idea onto a solid theoretical footing.
On the experimental side, the advent of the LHC has made it possible to produce bb abundantly, such that the production of e.g. the S-wave bound states in p + p as well as Pb + Pb has been determined with unprecedented precision [18]. While calibration of absolute yields is ongoing [19,20], the measured relative suppression of excited states in the presence of a deconfined medium already provides a first stringent benchmark for theoretical descriptions of in-medium bb.
The large mass of e.g. the ground state Υ also means that its formation time from the bb pairs created in the initial hard phase of a heavy-ion collision τ bb→Υ ∼ r Υ = 0.14fm [21] is much shorter than the QGP thermalization time. Of course the separation of initial production and QGP formation is not as pronounced for the excited states, such as χ b . Nevertheless we will proceed, based on the working hypothesis that fully formed bb bound states enter the quark-gluon plasma.
This scenario differs markedly from the description originally drafted for the case of charmonium, the bound states formed by a c quark and its anti-partner. Due to the smaller charm mass m c = 1.275 ± 0.025GeV (M S) [12] and correspondingly larger charmonium radii, it is expected that many cc pairs will be produced in the initial stages of the collision, which subsequently enter the QGP without having formed a fully bound state. In the presence of thermal fluctuations it is then argued that the strong interactions are not capable to support the formation of a bound ground state already at temperatures as low as T ∼ 1.2T C . Hence the suppression of the measured abundances was proposed as a prime signal of QGP formation.
This train of thought, first introduced in the seminal paper [22] and later refined e.g. in [23] relies on an investigation of the properties, such as the radii, of the Eigenstates of the in-medium Hamiltonian for the two-body cc system. These solutions describe quarkonium formed from cc pairs, which are essentially in thermal equilibrium with the surrounding medium, since they propagate in time simply through a phase factor.
A similar type of information is contained in the thermal spectral functions ρ(ω) calculated using the in-medium heavy quarkonium Schroedinger equation [24,25] or directly from lattice QCD [26][27][28]. In vacuum the presence of individual delta peaks clearly outlines possible bound states for the cc pairs to end up in. At finite temperature, broadening and disappearance of these structures is related to the fact that Debye screening and Landau damping [29] weakens the net attractive forces responsible to bind thermal qq pairs into a well defined heavy quarkonium.
To conclude as to whether the heavy quarkonium bound state will be able to form or not however relies on an inspection by eye. The absence of a clear concept of what constitutes a bound state in this otherwise theoretically solid assessment leads to ambiguities, reflected in different interpretations of the broadening, i.e. threshold enhancement vs. melting, discussed e.g. in [27,30].
In other words, while thermal spectral functions can help us to understand whether fully thermalized QQ pairs will be able to bind into quarkonium, we require different information to approach Bottomonium in a scenario assuming early-time formation and slow thermalization. The bb states then traverse the QGP truly as probes and not as part of the medium. Furthermore, as the bound states possess a long lifetime relative to the nuclear scales, their decay into gluons may occur at a point where the QGP has already ceased to exist. Both issues preclude us from making direct comparisons between the calculated Bottomonium spectral functions and the measured dilepton spectra by e.g. the CMS collaboration [18], which are related only if the dileptons are emitted from the annihilations of fully thermalized quark anti-quark pairs within a plasma [3,24].
Take the bb ground state in the vector channel, Upsilon, for example. When traveling in the QGP, it is not an Eigenstate of the in-medium Hamiltonian and thus will experience mixing with the excited states Υ and Υ , a phenomenon which has recently been emphasized in [1,[31][32][33][34]. In addition the thermal fluctuations in the medium will induce additional excitations and de-excitations between the different Bottomonium states along the time evolution. The language of open quantum systems and in particular the concept of stochastic potential will help us to describe this situation in an intuitive fashion.
In section II we briefly review the origin for the Schroedinger equation for heavy quarkonium and the inter-quark potential derived from the underlying field theory QCD via effective field theory methods. The fact that this potential takes on complex values at finite temperature is interpreted in terms of decoherence and we describe how the concept of stochastic potential allows us to set up an approximate description of the Bottomonium real-time dynamics. We emphasize that the thermal fluctuations and the imaginary part of the EFT potential are intimately related in this approach.
For a first and thus necessarily rather crude estimation on the evolution of Bottomonium in the QGP we carry out three dimensional simulations of the stochastic dynamics in section III. We assume a static medium described by the perturbative EFT potential [35] at T = 2.33T C and observe how individual states as well as non-thermal ensembles of Bottomonium evolve over time. We discuss the dependence of the abundances of heavy quarkonium states on information beyond what is contained in the thermal spectral functions and conclude in section IV with a list of necessary improvements to our method before a quantitative comparison to experiment can be attempted.
2 From complex to stochastic potential

A Schroedinger equation for heavy quarkonium
If there exists a separation of scales in a system and we are concerned with the physics only at one of these scales, the framework of effective field theory allows us to systematically focus on the energy range of interest. In the case of Bottomonium [14] we wish to understand the interaction of the medium gluons with the bb bound state which occur at the scale of the medium temperature T ∼ 200MeV, much smaller than the rest mass of the individual quarks m b = 4.66GeV [12]. Together with the fact the m b is much larger than the intrinsic scale of QCD Λ QCD ∼ 200MeV we expect that neither thermal nor quantum fluctuations will be able to pair create heavy bb pairs, so that a non-relativistic description is possible [13,36].
The starting point for the EFT [13,16,17] is a description of heavy quarkonium in the language of the fundamental field theory QCD, which relies on correlation functions of meson operators M (x, y, t) = Ψ(x, t)ΓU (x, y)Ψ(y, t). In anticipation of the test charge character of the heavy quarks, a finite separation is introduced by hand r = x − y, which requires us to connect the fields with a straight Wilson line U (x, y) if gauge invariance shall be preserved in the following. Γ denotes a Dirac matrix used to select the appropriate spin structure of the meson, e.g. Γ = γ µ for the vector channel particles.
The quantity we choose to help us understand the in-medium evolution of Bottomonium is the real-time forward propagator It encodes the correlations between a color singlet bare meson inserted into the medium at t = 0 and its counterpart after having evolved up to time t. While initial and final state are devoid of color, intermediate states contributing to Eq.(2.1) can have finite color charge, as medium gluons will be absorbed and reemitted over time. If D > (r, t) is evaluated as follows these medium interactions will be included.
To simplify the description we can exploit the separation of scales by integrating out from the above expression the hard energy scales of the rest masses. To this end we can expand the heavy quark action S ΨΨ in terms of increasing powers of m −1 q . To first order this leads to a decoupling of the upper and lower components of the original Dirac four spinor Ψ = (ξ, χ) into two two-component Pauli spinors representing the quark and anti-quark respectively. Rewritten in these non-relativistic degrees of freedom, we obtain the heavy quark action of the effective field theory NRQCD. Note that if the heavy mass expansion is continued to higher orders the residual influence of the higher lying scale eventually needs to be taken into account by correctly matching the Wilson coefficients that accompany each term of the NRQCD action [14].
The quadratic structure of the NRQCD action allows us to integrate out the heavy fermion fields explicitly [16,36]. Subsequently it becomes possible to leave the language of fields all together and represent the forward correlator in terms of quantum mechanical path integrals, where the fluctuating path z i (t) and the conjugate momenta p i (t) of the quark and antiquark are the dynamical degrees of freedom.
Note that for static quarks the term on the right is nothing but a Wilson loop in Minkowski time along the rectangular path the heavy quarks trace out in time. In that case the forward propagator obeys a Schroedinger equation and we can identify We find [37] that the potential term Φ(r, t) can actually be a time dependent quantity, which relaxes to a constant only at late times since the (instantaneous) effective potential V EFT (r) replaces a retarded interaction mediated by gluons.
Evaluating the real-time Wilson loop and thus V EFT (r) in perturbation theory [15,29,35], in AdS/CFT [38] and on the lattice [39][40][41] revealed that in general it takes on complex values at finite temperature. In the QGP, i.e. for T > T C the real-part shows the characteristic Debye screening behavior induced by the presence of deconfined color charges, as one would expect from the analogy with the electromagnetic plasma. On the other hand, since the medium partons are very light they can easily scatter with the gluons that mediate the strong force between the heavy quarks (Landau damping), which is reflected in the presence of an imaginary part that increases with temperature.
This tells us that the QCD derived potential is qualitatively different from purely real potential models, such as the color singlet free or internal energies [42][43][44], which have been used extensively before the effective field theory approach matured. I.e. even if the real part of the potential is quantitatively captured in a satisfying manner by the model potential [41,45], the absence of the imaginary part is a sign of missing physics.
The fact that the inter-quark potential defined from an effective field theory description is complex urges us to carefully examine its meaning. The central observation to make is that the Schroedinger equation (2.3) does not act on the wave-function ψ QQ (r, t) of the two-body system but instead operates on the correlation function D > (r, t). 1 Hence the presence of an imaginary part in V EFT (r) does not make any immediate statement on a reduced probability to find the heavy Quarkonium or in other words on the decay of the QQ pair. In particular the infinitely heavy quark and anti-quark pair for which the complex potential has been originally defined, will always remain present in the system, as its decay into gluons is fully suppressed in the non-relativistic description.
Im[V EFT (r)] = 0 instead reflects the fact that the thermal fluctuations of the medium decorrelate the current state from its initial conditions. Note that this too is not a statement on the melting of heavy quarkonium states, but only on the disappearance of correlations. To adequately take into account this phenomenon of thermal decoherence in the description of heavy quarkonium evolution, we turn to an approach based on an open-quantum systems viewpoint.

An Open-Quantum systems viewpoint
The idea of treating heavy quarkonium as an open-quantum system has received increasing attention in recent years [1,31,32,34,46,47]. The formalism originally developed in the context of condensed matter systems [48] is particularly suited to the study of Bottomonium since the intrinsic separation of scales allows a clear distinction between environment (QGP) and test system (Bottomonium) degrees of freedom. Our aim is to use the stochastic potential approach to heavy quarkonium evolution, first introduced in [1]. It makes close contact to the potentials derived in the EFT's and by doing so allows us to incorporate non-perturbatively information about the medium if e.g. potentials extracted from lattice QCD are used.
The overall system consisting of the environment, represented by the QGP, and the QQ test system is described by a hermitian Hamiltonian The states of the system evolve unitarily under H full and the density matrix of states σ(t) follows the von-Neumann equation. For simplicity we assume a separation of the initial states σ(0) = σ med ⊗ σ QQ (0). In the following we wish to track the evolution of the QQ subsystem by using only its own degrees of freedom , which is why we proceed by tracing out all the medium degrees of freedom. One arrives at the following heavy quarkonium density matrix where ψ QQ (r, t) denotes the wavefunction of heavy quarkonium two-body system in relative coordinates.
The phenomenon of decoherence formally tells us that [48] there exists an apriori unknown basis of states in which the reduced density matrix of states Eq.(2.7) becomes diagonal over time. The selection of this basis as well as the timescale for decorrelation depends on the interaction between the medium and the test-system. Intuitively decoherence it is related to the fact that thermal fluctuations distort the wavefunction of the test-system as it evolves in time, so that taking the thermal average over many realizations of the system leads to a loss of norm in the averaged wave-function.
We follow [1] in setting up an description of the time evolution of the static quarkonium system on the level of its wavefunction ψ QQ (r, t). Our goal is to both implement the concept of decoherence encoded in the imaginary part of V EFT (r) while still preserving the purely real values of the energy for the two-body system. In particular this prohibits us from directly using Im[V EFT ](r) in the Hamiltonian of the system. 2 A rigorous derivation of the stochastic potential, as well as its generalization to finite mass using the Feynman-Vernon influence functional has been worked out in ref. [47], .
The idea behind our approach is to reinterpret the presence of the imaginary part of the EFT potential V EFT (r) as an uncertainty in the values of a purely real potential V QQ (r) governing the evolution of ψ QQ (r, t). I.e. once a bound state of heavy quarks held together by a color string enters the QGP it will experience a weakening of the force between its constituents (Debye screening) and the gluons of the medium will continuously kick the color string (Landau damping) leading to slightly different values of the screened potential at subsequent times.
A naive way to formalize this idea is to construct a fully hermitian time evolution operator whose real potential term is disturbed by Gaussian, i.e. Markovian noise Θ(r, t) = 0. In order to encode the thermal properties of the QGP medium we allow the noise to carry a non-trivial spatial correlation structure Θ(r, t)Θ(r , t ) = 1 ∆t δ t,t Γ(r, r ). projec.on,of,, ψ QQ (r,t) , on,φ n (r) , ψ QQ (r,t),in,a,stochas.c,, Schrödinger$like,equa.on , Expanding Eq.(2.8) according to the rules of Ito calculus, we obtain the following stochastic Schroedinger-type equation Note that the term on the RHS with a factor i originates from the fluctuations Θ(r, t).
Evolving a single realization of the heavy quarkonium system, i.e. an individual wavefunction ψ QQ (r, t), along real-time still preserves its norm |ψ QQ |(t) = 1. Contrast this to the behavior after taking the thermal average is not yet enough to fully specify the dynamics of the system, indeed up to now, the off-diagonal parts of Γ(r, r ) are not determined. Since Eq.(2.10) does not contain further reference to the offdiagonal terms we might ask, whether such additional information is actually relevant to understanding Bottomonium evolution?
Up to now we have only spoken about the wavefunction of the heavy quarkonium system. To make a statement relevant for phenomenology, we however need to know the abundances of the different Bottomonium states present in the medium after a certain time t. In particular we ask (Fig. 1): what is the probability to find an eigenstate φ n (r) of the vacuum Hamiltonian H VAC QQ φ n (r) = E n φ n (r) at time t as it is the remnants of such states which are measured in experiment. In the stochastic potential approach the answer lies in the so called state admixture or survival probability of the state φ n (r) denoted by c nn (t). It is obtained by projecting the current wavefunction ψ QQ (r, t) onto the corresponding vacuum wavefunction φ n (r) c nn (t) = d 3 r d 3 r φ * n (r) Ψ QQ (r, t)Ψ * QQ (r , t) φ n (r ). (2.11) Note that this quantity depends on the evolution of the density matrix σ QQ (t, r, r ) = Ψ QQ (r, t)Ψ * QQ (r , t) whose time dependence is actually sensitive to the off-diagonal elements of Γ(r, r ) (see Sec. IV in [1]). This tells us that observing the evolution of Bottomonium gives us access to information that goes beyond what we can obtain from the EFT potential and thus by extension the thermal spectral functions.
The difference to inserting an imaginary part directly into the Schroedinger equation of ψ QQ (r, t) is now also evident. Instead of calculating the time evolution of the fluctuating ψ QQ (r, t) one actually determines the time evolution of ψ QQ (r, t) which however does not allow us to extract the time evolution of states in Eq.(2.11), since in general We have to caution the reader however that the intuitive setup presented here of course has its own deficiencies. The possibility to describe the heavy quarkonium evolution solely by a real stochastic potential has been shown to be valid only in the infinite mass case. By allowing a finite mass term in Eq.(2.8) we venture beyond its formal applicability, since the quantum analog of the classical drag force is not taken into account. The system described by Eq.(2.9) thus cannot reach thermal equilibrium even at asymptotically late times and will eventually accumulate energy at a linear rate in time [1,47]. Since our goal here lies in understanding relatively early time physics t ≤ 2fm where the artificial linear energy rise is not yet dominating, we believe that the numerical simulations presented in the following section do contain insightful information on the Bottomonium evolution.

Bottomonium melting via a stochastic potential
By constructing an approximate open-quantum systems description of heavy quarkonium based on the concept of stochastic potential in the previous section we are able to relate the real-and imaginary part of the EFT potential V EFT (r) to the heavy quarkonium potential V QQ (r) and the diagonal noise correlations of the medium Γ(r, r). In the following we will perform numerical simulations based on Eq.(2.9) and investigate the real-time evolution of Bottomonium at high temperatures to obtain a first view on the melting process from the stochastic potential approach.
We describe the vacuum properties of the Bottomonium system by utilizing a Cornell type potential in the Hamiltonian (2.9) as the wavefunction norm (red) remains at unit value and the energy of the system remains purely real. (blue) E[ψ] rises due to the influx of energy from the thermal medium but due to the absence of the correct momentum drag it will not settle to a time independent value and an artificial linear rise will dominate at later times [47]. The actual correlation function Γ(r, r ) obtained from sampling 60 noise realizations. Note that the diagonal follows the values of the imaginary part of the EFT potential Γ(r, r) = −Im[V EFT ](r) and the transversal falloff obeys a Gaussian Γ(r, r ) ∝ exp[−|r + r | 2 /l 2 ] with characteristic length l = 1/T . with the bare bottom mass m b , the string tension σ and the strong coupling α S chosen such that the ground state Upsilon has a mass of m Υ = 9.482GeV. Letting the linear term in Eq.(3.1) go over smoothly to a constant value at r SB = 1.32fm to take into account QCD string breaking leads to a spectrum in which there exist e.g. three bound S-states (Υ, Υ , Υ"), three bound P-states ( 1 χ b , 2 χ b , 3 χ b ) and two bound D-states ( 1 Υ 3 D 2 , 2 Υ 3 D 2 ) below the B-meson threshold.
For the description of the in-medium evolution we have to specify the real values of the heavy-quark potential V QQ (r) as well as the correlations of the thermal noise Γ(r, r ). Since in this exploratory study we will work at rather high temperatures T = 2.33T C = 396MeV we rely on the perturbative values of the EFT potential shown in Fig. 2. The functional form taken from [35] reads where we set the Debye mass to m D = 1GeV and the gauge coupling to g = 2.14.
To specify the off-diagonal elements of the noise we make the simple Ansatz that they decay in a Gaussian fashion with a characteristic correlation length given by l = 1/T (see e.g. Fig. 3) (3.5) While this choice seems reasonable at high temperatures, the behavior of the phenomenologically relevant QGP close to the phase transition (or rather cross-over) is characterized by growing correlation lengths and deviations from a simple Gaussian must be expected. Ultimately we will need to to find ways both conceptually and technically how Γ(r, r ) can be determined from an appropriate lattice QCD observable in its entirety. With these prerequisites set, we proceed by discretizing the system on a three dimensional hypercube with N = 256 points along each axis. The physical length of the box is chosen to be L = 7.68fm with a lattice spacing of a = 0.03fm so that all three bb S-wave bound state wavefunctions both fit into the available volume and are sufficiently well resolved (see gray curve on the left in Fig. 2). To numerically obtain these eigenfunctions for the vacuum potential of Eq.(3.1), which we will use as initial conditions for the in-medium evolution, we rely on the PETSC [51] and SLEPC [52] library.
The dynamics of the Bottomonium states according to Eq.(2.9) is subsequently determined using the unconditionally stable and norm preserving Crank-Nicholson method [53] with a time step of ∆t = 7.5 × 10 −4 fm. We prepare ten realizations of the system based on different initial conditions for the noise generator. In each run the norm of the wavefunction remains at unity as shown in the right panel of Fig. 2. For every fourth step in time we determine the admixtures of the individual bound states c nn (t) by projecting the stochastically evolved wavefunction onto the Eigenfunctions φ n (r) according to Eq.(2.11) and carrying out the noise average.

Pure initial state evolution
Before investigating the fully stochastic evolution of Bottomonium in a static QGP medium we first wish to observe its behavior in the absence of thermal noise. This case, which is known very well from the study of purely real model potentials (see also [1,[31][32][33][34]), already allows us to observe an important contribution to the dynamics: state mixing. The vacuum Υ initial state is not an Eigenstate of the Hamiltonian with Debye screened potential. Hence as shown on the left of Fig. 4, its admixture decreases steadily and instead the excited states  In both cases state mixing occurs since the vacuum initial state φ 0 (r) is not an Eigenstate of the in-medium Hamiltonian. Note that state mixing only changes the abundances of states with the same parity, i.e. no P,D etc. appears (left). Thermal noise contributes additional excitations and deexcitations to the time evolution which (albeit weakly) also excited states with different angular momentum.
Υ and Υ" become populated. At t ∼ 2fm almost equal amounts can be found. Note that there is a distinct time-lag between the increase in the numbers of Υ and Υ reflecting their different masses. Summing up the admixtures of the bound states we find that this value deviates slightly from one, the continuum states are only weakly populated. Note that because the Hamiltonian respects parity, we only find mixing between states that have the same angular momentum, i.e. no P or D states appear during the evolution.
In the case where we take the effects of thermal decoherence, i.e. of the imaginary part of the EFT potential into account, we still find that state mixing is an important part of the dynamics. As shown on the right of Fig. 4 the stochastic noise leads to a much faster depopulation of the Υ ground state. As the Θ(r, t) in our approach does not discriminate between states of different angular momentum it thermally excites all of the different P and D states albeit with a small magnitude, as well as the continuum states. Since we do not yet have a baseline for the overall production of bb in heavy-ion collisions we concentrate here only on the relative abundances between the three S-states plotted in the lower part of the figure.
The evolution of initial excited states plotted in Fig. 5 shows much shorter timescales for the initial depopulation as in the case of the ground state. As expected, if we start with Υ it affects both the bound states above and below, so that the admixtures of Υ and Υ" up to t < 0.8fm grow with similar speed. Since Υ" however cannot be efficiently fed by the ever declining Υ it subsequently depopulates again for t > 1fm. In the late time region t > 1.5fm where the quantitative reliability of our simulation is not guaranteed anymore, we find that the buildup in Υ will lead to a replenishment of the almost fully depopulated Υ state.
As last example in this subsection we take the evolution of an initial Υ" state which exhibits the most rapid depopulation of all three states. The total admixture of bound states exhibits two regimes of distinct rates of decrease in this case. The first is dominated by the faster than exponential decrease of the Υ , the second, slower one, follows the decay of the ground state which has taken over after t > 1.3fm.

Mixed initial state evolution
Since in the initial hard partonic stages of the collision at the LHC several bb pairs will be produced, we have to take into account the possibility that an ensemble of different Bottomonium states will enter the QGP together. Thus we have now two ensembles to take care of in the simulation, one concerns the realizations of the thermal noise, the other the probability to find a certain Bottomonium state as initial condition.
Here we make a naive ansatz and roughly take the relative abundances of the Bottomonium S-wave states obtained from the measurements carried out by the CMS collaboration [18] for p + p collisions at √ s = 2.67TeV. Using the ratio Υ : Υ : Υ = 5 : 3 : 2, we investigate the evolution of the admixtures over time, as shown in Fig. 6. The first observation is that at small times t < 1fm the excited S-wave states feed the ground state and lead to its relative enhancement. This timescale aligns reasonably well with the rapid depopulation of the Υ and Υ states observed in Fig. 5. If we take a look at the evolution at later times, a replenishment of the excited states occurs with Υ being the main benefactor as it is directly sourced by the ground state.  Figure 6. The evolution pattern for the Bottomonium family in the presence of a mixed initial state. Here we use relative abundances inferred from the CMS dilepton spectrum in p + p, i.e. Υ : Υ : Υ" = 5 : 3 : 2 as starting point for the evolution in a medium at T = 2.33T C . In this highly idealized setting, we observe that the interplay between replenishment and suppression of the individual Bottomonium states for times t < 1.5fm actually leads to a slight enhancement of the ground state compared to the excited states.
While it would be rather bold to claim any quantitative relation of this naive calculation with measurements in experiment it is reassuring that even in this highly idealized setting one observes at time t < 1fm a behavior akin to melting of the higher excited states, while populating the ground state.

Discussion
The results for the evolution of a mixed initial state are encouraging in that they hint at excited state melting, obtained from a consistent potential-based description of heavy quarkonium including the effects of thermal decoherence while circumventing the interpretational difficulties associated with complex system energies. Nevertheless there is still a long way ahead of us before we can quantitatively compare to experiment or as would be the final goal to infer novel information about the medium from the measured Bottomonium di-lepton spectra.
The first reason is that we do not yet have any reliable knowledge about the offdiagonal elements of the noise correlations. All the results shown above were based on a phenomenologically motivated choice of correlations, decaying according to a Gaussian with correlation length l = 1/T . At high temperatures with short correlation lengths such behavior might be plausible, close to the phase transition it is far from certain that it persists. Sadly we cannot ignore this issue as changing the values of l, while not reflected on the level of the average wave-function or equivalently the thermal spectra, leads to discernible changes in the state evolution shown in Fig. 7. It would be very helpful to find a field theoretical object which can serve as its counterpart in an EFT description. Just as the forward correlator allowed us to extract V QQ (r) and Γ(r, r), an appropriate four point meson correlator might supplement our knowledge on Γ(r, r ) for r = r . The second and more fundamental reason is the inability for Eq.(2.9) to correctly describe the thermalization of the heavy quarkonium with its surroundings. As was shown by Akamatsu in [47], if we wish to extend the predictive range of our computation to late times t ∼ 10fm we will have to make explicit the dissipative character of the dynamics by distinguishing forward and backward propagating modes and simulate them individually.
The conceptually improved description [47] takes into account the color degrees of freedom explicitly, reflected e.g. by the presence of colored noise, so that both the propagation of individual quarks and their bound states can be captured in medium. I.e. both color singlet and octet configurations are dynamical degrees of freedom. Inspection of their equation of motion hints at another limitation of the simple color neutral stochastic potential used here. As long as the heavy quarks in a singlet configuration remain close enough together, the gluons will not excite them into an octet state and we can safely assume as we did that the singlet propagates into a singlet. If however they become separated such that they appear as individual quarks to the medium, i.e. further than the thermal correlation length, it will quickly color rotate them into an octet configuration which then would have to be included explicitly.
The choice of initial conditions might also be ameliorated using the recently presented precision measurements for the χ b states [54] in p + p collisions. It will be interesting to see how these states with finite angular momentum evolve in the medium and whether their presence in a mixed ensemble does affect the S-wave abundances through (de)excitations induced by the stochastic noise.
Future simulations should also consider that the QGP actually evolves rapidly. It expands and cools, as can be explored by relativistic hydrodynamics simulations. While initially restricted to perfect liquids, fully causal approaches to viscous fluids are nowadays  Figure 8. (left) Bottomonium evolution from an pure initial Υ state based on an explicit imaginary part of V QQ (r) in the Schroedinger equation (2.9) without stochastic noise Γ(r, r ) = 0 . As in the case for a purely real Hamiltonian we do not find any excitations of states with angular momentum different from zero. Comparing to the stochastic evolution in the right panel of Fig. 4, we find that e.g. the melting of the ground state appears to proceed more efficiently here. The damping of the wavefunction will ultimately lead to vanishing admixtures at late times and not to a thermal distribution. (right) Decay of the wave function norm over time (red solid) and the accompanying decrease in the values of the system energy (blue). Note that a finite but small imaginary part (magenta) is induced in the energies of the system and that the energy monotonously decreases. This reduction of the real-part of the energy does not comply with intuition. If a vacuum state is inserted into a thermal medium we would expect that the kicks of the medium constituents will actually transfer energy into the system until the a common temperature is reached eventually and the energy becomes time independent.
available [55,56]. On the other hand progress has been made by constructing analytic solutions to the hydrodynamic equations [57] that can take into account an equation of state provided by lattice QCD. After checking that the potential based approach remains valid in the presence of a rapidly changing temperature, the stochastic simulations should be amended by changing the parameters of the potential accordingly.
Another interesting feature to study within the stochastic approach is the momentum dependence of the heavy quarkonium melting process. At early times, i.e. before momenta are able to equilibrate, the QQ state is certainly not at rest and it traverses the QGP. Studies of thermal spectral functions in lattice NRQCD [58] as well as using EFT methods [59] revealed that its probability to remain in a well defined bound states is indeed affected by the velocity of the constituents. Since the potential based approach allows to equip the vacuum wave function with a finite momentum by a simple phase factor we should aim at simulating Bottomonium moving through an appropriately cooling medium.

Summary and Conclusion
After revisiting the derivation of the Schroedinger equation (2.3) and the associated effective field theory potential (2.5) for the forward correlator of heavy quarkonium, we elaborated on the presence of an imaginary part of V EFT (r) in the QGP in terms of thermal decoherence.
In order to consistently incorporate the effects of decoherence in a potential based, i.e. nonrelativistic but otherwise dynamic, real-time description of Bottomonium wavefunctions (not correlators), we followed the Open-Quantum systems viewpoint introduced in [1].
Constructing an approximate unitary time evolution operator (2.8) for the heavy quarkonium wavefunction ψ QQ (r, t), based on a real but stochastic potential, we can calculate the admixtures c nn (t) of individual vacuum states φ n (r) as the system propagates in medium by appropriate projections (2.11). This approach invites the identification of the imaginary part of the EFT potential with the diagonal noise structure of the medium induced stochastic noise Γ(r, r), offering a (possibly non-perturbative) bridge between the EFT and the open-quantum systems language [1]. It differs thus from previous potential based studies (e.g. [33,49,50]) in that no explicit non-hermiticity enters the heavy quarkonium Hamiltonian, all Energies remain real-valued and actually increase over time (Fig. 2 vs. Fig. 8) as expected from inserting a vacuum state into the medium.
The dependence of c nn (t) on the density matrix of states and thus on the off-diagonal elements of the noise correlations (see Fig. 7) entails that by observing Bottomonium states over time we can access information beyond what is encoded in the EFT potential, and hence the thermal spectral functions, which can be determined from V EF T (r). I.e. if we already knew the temperature of the medium created in relativistic heavy ion collisions and a reliable EFT potential is used, the measured abundances of Bottomonium might contribute to shedding light on the correlation structure of thermal fluctuations within the QGP.
As a first attempt at estimating the behavior of different Bottomonium vacuum states as they propagate through the QGP, we investigate a highly idealized static thermal medium. At a temperatures of T = 2.33T C the perturbative EFT potential by Laine is used to set the parameters V QQ (r) and Γ(r, r), while for the off-diagonal entries we make a simple Gaussian ansatz with a characteristic correlation length of l = 1/T ( Fig. 2 and Fig. 3).
As was known from the study of model potentials and has been emphasized recently in [1,[31][32][33][34], we find that the mixing of different vacuum Eigenstates due to the presence of a screened real part of the heavy quark potential V QQ (r) is a significant contribution to the overall dynamics (Fig. 4). Note that in the absence of thermal noise (i.e. also for purely real EFT potentials), only states with the same angular momentum as the initial state participate, which however also includes unbound continuum states.
In the presence of thermal fluctuations, i.e. when taking into account the imaginary part of the EFT potential, the admixture of bound states decreases much more rapidly. Starting from pure initial states ( Fig. 4 and Fig. 5), one finds that an intricate interplay between mixing and the stochastic (de)excitation of states leads to non-trivial patterns of decay and replenishment of the S-wave states. Since the noise can also excite states of different angular momentum, we observe a finite but relatively small admixture of P, D, ... states. Until the open-heavy flavor measurements currently performed at the LHC have established a baseline for the overall number of bb produced in Pb + Pb collisions we refrain from absolute comparisons and focus on the relative abundances of the individual states.
Starting from a mixed initial state consisting of Υ : Υ : Υ" = 5 : 3 : 2 ratio of S-wave states, inspired by the dilepton spectra determined by the CMS collaboration in p + p collisions, a mild overall melting effect of the excited states is obtained at early times t < 1fm (Fig. 6). We find a tendency for regeneration of the excited states at later t > 1.5fm which is to be expected from the behavior of the pure initial state runs. Due to the very naive setting of Bottomonium at rest in a static thermal medium this result should not be thought of as a quantitative statement but rather as a motivation to develop further the Open-quantum systems approach to heavy quarkonium in general and the stochastic potential framework in particular, aiming at more realistic scenarios.
The first step will be to include the effects of an evolving temperature obtained from numerical or analytic solutions of the hydrodynamic equations. Together with endowing the initial wavefunction with a finite momentum and taking into account more than just the initial S-wave admixtures it should be possible to obtain a more realistic estimation of the Bottomonium evolution at early times. The implementation of the improved formulation of the stochastic dynamics based on the Feynman-Vernon influence functional should concurrently be pursued to enable a more reliable investigation of the late time approach to thermal equilibrium which we expect will be of use in building a bridge from the language of thermal lattice QCD spectral functions to the phenomenology of Bottomonium melting.