Classical model for diffusion and thermalization of heavy quarks in a hot medium: memory and out-of-equilibrium effects

The description of high energy collisions as the interaction of sheets of color-glass-condensate (CGC) [1–8] is one of the most intriguing approaches to the study of relativistic proton-proton (pp), proton-nucleus (pA) and nucleus-nucleus (AA) collisions at the Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC) energies. One of the predictions of the CGC effective theory is the formation of strong gluon fields in the forward light cone named glasma [9–20], which consists of the out-of-equilibrium longitudinal color-electric and color-magnetic fields. These are characterized by large gluon occupation number, $ A_\mu^a \simeq 1/g $, with g representing the strong coupling, so that they can be described by a classical field theory, namely the classical Yang-Mills (CYM) theory. Heavy quarks are excellent probes of the system created in high energy nuclear collisions, both of the pre-equilibrium phase and of the thermalized quark-gluon plasma (QGP), see [21–49] and references therein. Their formation time is very small compared to the light quarks, $ \tau_{\rm form} \leqslant 0.1 $ fm/c for both charm and beauty quarks. Since heavy quarks are produced immediately after the collision, they can propagate in the evolving glasma fields and probe its evolution. Propagation of heavy quarks in glasma has been studied previously with the Fokker-Planck equation [47], as well as the Wong equations coupled to the evolving glasma background field [48, 49]. The common result of these references was the finding of the diffusion of heavy quarks in the dense gluon medium, very similar to that in a thermal medium. On the other hand, a substantial drift was not observed in [48, 49] in the kinematic regime examined, which paved the way for a possible solution of the $ R_{AA}/v_2 $ puzzle of heavy quarks in Pb-Pb collisions [43].

The aim of the present work is to use a simple classical model to describe the motion of heavy quarks in a dense and hot bath, that has been used extensively in connection to the dissipation in quantum Brownian motion, see for example [50–53] as well as references therein. The advantage of this model, despite its simplicity, is that it has the potential to describe several aspects that might be important for understanding the drift and diffusion of the quarks, among them the propagation in a bath out-of-equilibrium like glasma, as well as the memory effects in the dissipation kernel. In the particular context of high energy nuclear collisions, memory has started to attract some interest [44, 45, 54–56] , as it might affect the di-lepton yields [57] as well as radiative energy loss [58]. In this model, the bath is described as a system of harmonic oscillators. The initial conditions for the oscillators are randomly assigned. The expectation values for the squared momenta and positions are in agreement with the equipartition theorem in the case of a bath in thermal equilibrium, or to the energy per particle in the case that the bath is out-of-equilibrium. To model an out-of-equilibrium bath, we use a distribution inspired by the CGC models, which is flat below a certain energy scale, the saturation scale $ Q_s $, and then rapidly drops to zero [59–61]. The coupling of the bath to the heavy quark is pretty simple, $ H_{\rm int}\propto Q \sum_k x_k $, where Q is the generalized coordinate of the heavy quark and $ x_k $ denotes the position of the $ k^{\rm th} $ oscillator. This simple coupling is sufficient to emphasize the fact that in order to achieve energy loss during Brownian motion, the back-reaction of the oscillators on the bath needs to be taken into account. In fact, a heavy quark produces a disturbance that is superimposed on the harmonic motion of the oscillators, and this disturbance couples to the heavy quark itself causing dissipation.

One of the ingredients of this model is the distribution of oscillators, $ {\rm d}N/{\rm d}\omega $, where $ \omega $ denotes the proper frequency of the oscillator. This quantity affects the shape of the dissipative kernel. In particular, if $ {\rm d}N/{\rm d}\omega\propto\omega^2 $ then the kernel is of the Markov type, $ \gamma(t-t^{\prime})\propto \gamma\delta(t-t^{\prime}) $ , and there is no memory. On the other hand, by deforming the above distribution it is possible to obtain a kernel with memory, making it possible to study the effect of the latter on physical processes like thermalization. This model allows to obtain almost all analytical results, albeit in the weak coupling regime, allowing to understand how the distribution of the bath and the memory kernel affect a physical process.

The main difference between our study and the studies based on generalized Langevin equations is that in our study the motion of the oscillators and of the heavy quark is deterministic, while in the Langevin equations the noise is a random number with a given correlation. In the present model, the random nature of the bath comes into the game as soon as we assume a distribution for the initial conditions of the oscillators. Averaging the physical quantities over the initial conditions amounts to taking the ensemble average for the system. The apparent merit of the present model is that it is very easy to implement a bath that is out-of-equilibrium, which is an important aspect to consider when studying the diffusion of heavy quarks in the evolving glasma. Moreover, the dissipative kernel appears clearly as an effect of the coupling of the heavy quark to the disturbance it creates in the bath, and this might be useful for implementation of a dissipative term in field theories.

The plan of the article is as follows: in Section 2, we describe the model in detail and derive the evolution equation for the heavy quark momentum, P, and define the shape of the dissipative kernel. In Section 3, we solve exactly the evolution equation for P and estimate the thermalization time, and study how memory affects this quantity in the model. In Section 4, we briefly describe the equilibrium state. Finally, in Section 5 we present our summary and conclusions.

We consider the classical dynamics of a nonrelativistic particle with mass M interacting with a bath of $ N\gg1 $ harmonic oscillators characterized by mass m and frequencies $ \omega_k $. We also introduce a distribution function for the bath of oscillators, $ {\rm d}N/{\rm d}\omega $, with a normalization

In the above equation we have assumed that the oscillator density is an even function of $ \omega $. As will be clear soon, $ {\rm d}N/{\rm d}\omega $ (or to be precise, its Fourier transform) is the only quantity needed to build the dissipative kernel. Changing its analytical form allows to switch from a kernel without memory to one with memory. The Hamiltonian is [50]

where $ M\gg m $, and Q corresponds to the coordinate of the heavy quark. The interaction depends on the coupling g, which has the dimensions of a cubic energy, that we take to be the same for all oscillators. In order to mimic a hot QCD medium it would be desirable to have a coupling that depends on the energy scale of the bath, namely on the temperature. While its implementation would be straightforward, we leave it to a future work because the main calculations of equilibration are unaffected by this additional complication. The term $ {\cal H}_R $ in Eq. (2) is a renormalization potential [53] that is introduced to balance the effect of the back-reaction on the bath, leaving $ V(Q) $ as the potential. We will specify this term later and its purpose will appear less obscure as soon as the dynamical equation for P is derived. In actual calculations we replace summations over oscillators with an integral: for any f we use

For a bath in thermal equilibrium, the initial conditions for the oscillators are distributed randomly according to the standard Boltzmann distribution at temperature T,

where we have used

which imply $ \langle p_{0k}^2\rangle = mT $ and $ m\omega_k^2\langle x_{0k}^2\rangle = T $, in agreement with the equipartition theorem. On the other hand, when considering an out-of-equilibrium system we introduce a Fermi-Dirac-like distribution, namely

where $ {\cal T} $ and $ {\cal M} $ are two parameters with the dimension of energy. When $ {\cal M}/{\cal T}\gg1 $, the out-of-equilibrium distribution becomes a theta-function that has been used in literature to mimic the gluon distribution in glasma [59–61], with $ {\cal M} $ playing the role of the saturation scale $ Q_s $. For this reason, we call Eq. (7) the CGC-like distribution. The form of Eq. (7) also keeps $ x_{0k} $ and $ p_{0k} $ uncorrelated. While it would be interesting to study the effects of correlation between coordinate and momentum space, we neglect them for simplicity. The bath distribution in Eq. (7) implies $ \langle p_{0k}^2\rangle = m{\cal T}_{\rm eff} $ and $ m\omega_k^2\langle x_{0k}^2\rangle = {\cal T}_{\rm eff} $ , where we have defined the effective temperature of the bath,

where$ z = {\rm e}^{{\cal M}/{\cal T}} $ and $ {\rm Li}_s $ denotes the Polylogarithm function of order s. For example, in the case of a flat distribution $ {\cal M}/T\rightarrow \infty $ we find

Another ingredient of the present model is the distribution of the frequency of the oscillators, $ \omega_k $. We distribute it with the density $ {\rm d}N/{\rm d}\omega $ , which will be specified later, case by case, and which is normalized according to the condition in Eq. (1). We anticipate that the memory effect of the dissipative kernel is determined solely by the form of $ {\rm d}N/{\rm d}\omega $. Although this should be obvious, we remark that we have made an assumption of a rigid $ {\rm d}N/{\rm d}\omega $, meaning that although we consider the back-reaction of the heavy quark on the bath of oscillators, our dynamics is such that this back-reaction affects only the motion of the oscillators and not their internal structure, namely $ {\rm d}N/{\rm d}\omega $ is unaffected by the presence of the heavy quark.

The derivation of the dynamical equation for P is straightforward and well-known in literature. However, in order to emphasize that dissipation comes solely from the back-reaction of the oscillators on the bath, we review the main steps of the derivation below. This will also help to clarify the renormalization of the potential in Eq. (2). The starting point is the equation of motion of the oscillator k, that can be obtained easily from the Hamilton's equations, that is

This equation can be easily solved exactly using the retarded Green's function, namely

where we have defined

with $ {\cal G}_R^k(x) $ denoting the retarded Green's function of the oscillator,

Moreover, we have defined

which corresponds to the unperturbed harmonic motion of the oscillator (namely to the solution of the homogeneous equation), where $ \alpha_k = g x_{0k} $, $ m\omega_k\beta_k = gp_{0k} $ and $ x_{0k} $ and $ p_{0k} $ denote the initial conditions for coordinate and momentum of each oscillator in the bath. More generally, $ x_k^h(t) $ corresponds to the solution of the homogeneous equation of the oscillator. In this work we assume a statistical distribution of the initial conditions, but the dynamics of the heavy quark, as well as that of the oscillators, is purely deterministic. This approach is very similar to what is done within the CGC effective theory, in which a random distribution (usually Gaussian) of color charges is assumed, and the initial color fields are then computed and the evolution is studied by solving the CYM equations, which are deterministic. Within the CGC formalism the ensemble average is taken numerically by averaging the physical quantities over many different initializations of the color charges. In our work we follow a similar path. We first distribute the initial conditions of the oscillators by virtue of a statistical distribution, either at equilibrium or out-of-equilibrium. The evolution of the system is then studied with the classical equations of motion, and the ensemble average is obtained by averaging the classical solutions with the statistical distribution of the initial conditions. We remark that Eq. (11) is exact and contains the back-reaction induced by the interaction of the heavy quark with the oscillator, the latter being represented by the convolution of the Green's function with the generalized coordinate Q. The analytical form of the Green function will be specified later.

The equation of motion for P is just the second Hamilton's equation that can be read from Eq. (2), namely

Defining the source term, $ \xi(t) $, as

where $ x_k^h(t) $ denotes the harmonic motion of the momentum of the $ k^{th} $ oscillator, and taking into account Eq. (11) we can write

The convolution term on the left hand side of the above equation represents the coupling of the heavy quark to the perturbation of the harmonic motions of the oscillators, which is induced by the heavy quark itself. It is therefore the back-reaction term. Note that by integrating by parts we can write this convolution as

where $ F(x) $ is the primitive of $ G_R(x) $ , and we have taken into account that $ {\rm d}G_R(x)/{\rm d}x = -{\rm d} G_R(t-t^{\prime})/{\rm d}t^{\prime} $. The full equation for P can thus be written as

The role of $ {\cal H}_R $ appears now clearly: the term $ -F(0) Q(t) $ in Eq. (19) corresponds to a harmonic shift of $ V(Q) $ induced by the back-reaction. Taking

amounts to canceling this shift when the summation over the oscillators is performed, and $ V(Q) $ remains as the only potential. Using Eq. (20), we get

We now make the assumption that $ Q(0) = 0 $. This assumption can be relaxed easily, as will be clear when we present the exact solution of the equation of motion. Moreover, here we want to study solely the motion of the heavy quark in the bath, therefore we put $ V(Q) $ = 0. In this case we get

This equation looks like the generalized Langevin equation. However, we remark that the term $ \xi(t) $ on the right hand side of the above equation is not noise but a deterministic source, therefore the motion of the heavy quark is completely deterministic. An ensemble average is taken of $ P(t) $ , averaging over the initial conditions of the harmonic oscillators. We will prove later that for a Markov process the time correlator of $ \xi(t) $, defined in terms of the ensemble average, is indeed proportional to the $ \delta- $function in analogy with the Gaussian noise used in standard Langevin equations.

The last step is to write explicitly the function $ F(x) $. To this end we need the Green's function of the harmonic oscillator,

the primitive of $ G_R^k $ is easily calculated. Defining the dissipative kernel as

we can write the dynamical equation for P as

The derivation of Eq. (25) shows that the effect of coupling of the heavy quark to the oscillators appears in two places: on the right hand side of the equation, describing the coupling to the unperturbed harmonic motion of the oscillators, which is responsible for momentum diffusion as we will see shortly, and in the dissipative term on the left hand side, describing the coupling of the heavy quark to the perturbation created by its motion in the bath, which is responsible for dissipation. Neglecting this back-reaction amounts to neglecting the energy loss, and thermalization of the heavy quark will not be achieved. The back-reaction term appears as the next-to-leading order interaction. However, the drag and diffusion coefficients are of the same order in coupling, as we will see in the next subsection. Nevertheless, we will also provide an estimate of the kinematic regime in which neglecting the dissipative term is a reasonable approximation. In the continuum limit, using the replacement (3), we can write Eq. (24) as

which gives the Markov limit $ \gamma(t-t^{\prime}) = 2\gamma\delta(t-t^{\prime}) $ for $ {\rm d}N/{\rm d}\omega \propto \omega^2 $. In this case we get

We note from Eq. (26) that, apart a multiplicative constant, $ (1/\omega^2){\rm d}N/{\rm d}\omega $ is nothing else but the Fourier transform of the dissipative kernel. Therefore, it appears clearly that the shape of the latter can be modified by the former. The ensemble average is taken of the solution of Eq. (25). Practically this means that we are solving the deterministic equation of motion (25) for a set of initial conditions of the oscillators, and then averaging by means of the distributions specified in the previous subsection.

We now show that as long as we use Eqs. (4) and (7) to distribute the bath of harmonic oscillators, the dissipative kernel and the correlator $ \langle\xi(t)\xi(t^{\prime})\rangle $ satisfy the fluctuation-dissipation theorem regardless of the particular $ {\rm d}N/{\rm d}\omega $ chosen. To this end, taking into account that from Eqs. (4) and (7) it follows that $ \langle\alpha_k \alpha_{q}\rangle\propto \delta_{kq} $, $ \langle\beta_k \beta_{q}\rangle\propto \delta_{kq} $, $ \langle\alpha_k \beta_{q}\rangle = 0 $, we get

The above equation is general. If we now make the further assumption that

we get

Assuming that $ \langle p_{0k}^2\rangle $ is independent of the oscillators, which happens for the distributions used in this work, we can write

where now $ \langle p_{0}^2\rangle $ corresponds to the common average of the initial squared momentum of the oscillators.

A comparison of Eq. (31) and Eq. (24) shows that in this model the fluctuation-dissipation theorem is satisfied in the form

This relation is valid regardless of the particular form of $ {\rm d}N/{\rm d}\omega $, as long as the condition (29) is satisfied.

We remark that in order to satisfy the fluctuation-dissipation theorem the assumption (29) is crucial. If this condition is violated then the second term on the right hand side of Eq. (28) does not necessarily cancel out, and this might lead to the violation of the theorem. This might happen, for example, if in the initial conditions the degrees of freedom of the oscillators are thermalized at different temperatures. This situation represents a different out-of-equilibrium condition, and in the presence of interactions among the modes of the oscillators this condition should be removed by the dynamical evolution. We leave this possibility for a future investigation. In this study the distributions of the oscillators always satisfy the condition (29). Therefore, the aforementioned term cancels out and the fluctuation-dissipation theorem is always satisfied.

We now show that in this model the heavy quark experiences pure diffusion of the initial value of momentum in the case the dissipative term is neglected in a Markov process. Indeed, in this case the equation of motion is

Squaring and taking the ensemble average leads to

where we use $ {\cal E} = T $ for the system in thermal equilibrium, and $ {\cal E} = {\cal T}_{\rm eff} $ for the CGC-like distribution. The time dependence of the momentum spread comes from the integral over the frequencies of the oscillators on the right hand side of the above equation. In the case of a Markov process this is easy to compute since in this case $ {\rm d}N/{\rm d}\omega\propto\omega^2 $ , which implies

where we have defined the diffusion coefficient of kinetic energy as

The above equations describe standard diffusion with linear momentum spread with time.

We can compare the diffusion coefficient with the drag in the Markov limit. Indeed, from Eq. (26) we get

and comparing this with Eq. (37), we find

which corresponds to the fluctuation-dissipation theorem in the case of a Markov process.

In this study we consider a dissipative kernel with memory. As specified above, a Markov kernel (namely, a kernel without memory) corresponds to the limit $ \gamma(t-t^{\prime}) = 2\gamma\delta(t-t^{\prime}) $ . According to Eq. (26), this can be implemented by means of $ {\rm d}N/{\rm d}\omega \propto \omega^2 $. In order to keep the expressions that are easy to manipulate we consider a simple Gaussian form

which is normalized according to Eq. (1). In this equation we introduce the parameter a, which has the dimension of energy, that regulates the shape of the kernel in the frequency space. Using Eq. (40), we get from Eq. (26)

where

and

This form has the advantage of easily identifying the part of the kernel that gives the $ \delta- $function in the limit $ a\rightarrow\infty $ , since $ \Phi_a(x)\rightarrow\delta(x) $ in this limit. The form in Eq. (41) also suggests that the proper way to take the limit $ a\rightarrow\infty $ is to keep $ g^2 N/a^3 $ fixed, so that the overall constant in Eq. (41) is unchanged as the dissipative kernel is continuously deformed from the one with memory to the Markov one, $ \gamma(s) = 2\gamma\delta(s) $. In addition, we note that Eq. (41) suggests that we can define an effective coupling, which contains the information about the bath of oscillators and is independent of the mass of the heavy quark. As a consequence, once we fix M we can study the effect of a continuous deformation of the dissipative kernel by keeping $ \alpha^2 $ fixed. This limiting procedure gives meaningful results if we take both limits $ N\rightarrow\infty $ and $ a\rightarrow\infty $ .

The dynamical equation for P, Eq. (25), is an integro-differential equation that can be solved analytically by means of the Laplace transform. We present this solution first, and then discuss how it can be used to define the thermalization time, which we estimate for a given dissipative kernel.

In order to solve Eq. (25) we introduce the Laplace transform, namely

It is then straightforward to write Eq. (25) as

Here, $ \Xi(s) $ and $ \Gamma(s) $ correspond to the Laplace transforms of the source term and of the dissipative kernel, respectively. From Eq. (45), we can write the solution in the time domain as the inverse Laplace transform of $ \tilde{P}(s) $, that is

where $ \sigma $ is a real constant that is larger that the real part of all singularities of the integrand. Equation (46) corresponds to the exact solution of Eq. (25). In principle, once the dissipative kernel and the source are given, one can obtain the evolution of P by performing the integral above. The ensemble average can then be taken by averaging over the initial conditions of the oscillators in the bath.

Before applying the solution (46) to some specific dissipative kernel and source, we note that we can split this solution into two parts:

In this equation, $ P(0) $ is the initial value of the heavy quark momentum. The first term on the right hand side is therefore the only part of the solution that contains information on the initialization. Depending on the distribution of the poles of the integrand, the magnitude of this term can grow or be damped exponentially, as well as oscillate in time. If all poles have negative and nonzero real parts then the first term decays exponentially, implying that after some time the heavy quark loses memory of its initial condition. Its properties at asymptotic time depend only on the bath in which the quark propagates. This offers a criterion for thermalization, or more generally for equilibration, of the heavy quark with the medium. Indeed, we can define the thermalization time as the time needed for the quark to forget about its initial condition. This time, that we denote by $ \tau_{\rm eq} $, can be computed by studying the distribution of the zeros of $ s+\Gamma(s) $ in the complex plane Assuming the poles $ s = -\kappa_r + i\kappa_i $, with $ \kappa_r $ and $ \kappa_i $ corresponding to the real and imaginary parts of s, respectively, and assuming that all $ \kappa_r $ are positive, we can identify $ \tau_{\rm eq} $ with the smallest $ \kappa_r $. We note that using this definition, $ \tau_{\rm eq} $ depends only on the shape of the dissipative kernel in the frequency space and not on the initial conditions of the oscillator bath. Given a density $ {\rm d}N/{\rm d}\omega $, the equilibration time only depends on the density as well as on the coupling and the masses of the heavy quark and of the oscillators, but not on the distribution of the oscillators in the bath.

We now apply the solution (46) to the case of the Gaussian dissipative kernel Eq. (41). After writing down the exact solution, we use it to determine the equilibration time, as defined in the previous subsection, as a function of the shape parameter of the kernel. The Laplace transform of the kernel is

where N is the number of oscillators in the bath and we have used

We recall that we obtain the Markov limit when $ a\rightarrow\infty $ in the above equation. The Laplace transform of Eq. (16) is given by

Therefore, we can write the exact solution in the form

A quick calculation of the residues at infinity of the integrands in Eq. (51) at t = 0, which correspond to the integrals at t = 0, shows that the solution (51) is consistent with $ P(t = 0) = P(0) $: the second term on the right hand side of the equation vanishes at t = 0, while the first one is simply equal to $ P(0) $.

As discussed in the previous subsection, in order to estimate the equilibration time we have to study the distribution of the poles of the integrands on the right hand side of Eq. (47). The equation $ s+\Gamma(s) = 0 $ , with $ \Gamma(s) $ in Eq. (48), admits both negative real, $ s^\ell = -\chi^\ell $, and a set of infinite conjugate complex solutions $ s^n = -\kappa^n_r\pm i\kappa^n_i $ with $ \kappa^n_r>0 $. Besides, there are conjugate poles $ s_k = \pm i\omega_k $ in the second integral. We were able to find explicit solutions for the equation $ s+\Gamma(s) = 0 $ analytically only for small $ \alpha^2 $, in the limit of large a , that is for a bath that is not very far from the Markov limit. They were found by using the Newton-Raphson method iterated up to the second order (we checked that increasing the order of the iteration does not change the result in the weak coupling limit). These solutions are given by

with $ \kappa_i^n = -na $. Clearly the smallest zero is $ \chi^1 $ , and in the limit of large a, $ \kappa_r^n $ is very large, which means that the contribution of these poles is suppressed as soon as $ 1/\kappa_r^1 < t $ . Hence, in this time range we can limit the solution only to $ \chi^1 $. In this limit, by a straightforward application of the residue theorem we get

where c.c. denotes complex conjugation. The last line in the above equation corresponds to the contribution of the poles $ s = \pm i\omega_k $ and represents the heavy quark momentum at a very large time, while the first two terms are important only in the transient region.

As anticipated, part of the solution (54) depends on the initial conditions (the first line on the right hand side of the above equation). We identify $ \chi^1 $ with the inverse of the equilibration time, getting

The above equation shows that the in this model, the effect of memory in the dissipative kernel is to lower the equilibration time of the heavy quark when the effective coupling $ \alpha^2 $ is kept fixed. We were unable to find a detailed explanation of why this happens. A naive interpretation is that adding memory to the dissipative kernel amounts to lowering of the value of a and, because of Eqs. (40) and (43), keeping $ \alpha^2 $ fixed implies that the density of oscillators at a given frequency $ \omega $ is larger, resulting in a more efficient interaction with the medium and in a faster thermalization. To this end we remark that we have obtained Eq. (55) using a Gaussian kernel. It would be interesting to investigate other analytical forms of the dissipative kernel in order to verify if the lowering of the equilibration time is a more general feature, or if it is related to the specific form of the kernel used in our work. Furthermore, it would be interesting to study the impact of the memory effect on heavy quark observables at the RHIC and LHC energies, as well as on the heavy quark thermalization time, which is neglected in several recent calculations [37–42, 46]. It should be mentioned that the memory effect has an important influence on the final dilepton yields [57] and radiative energy loss [58]. This indicates that an analysis which describes the data correctly demands a proper inclusion of the memory effects (for more information see also [44, 45, 54–56]).

We note that the explicit distribution of the bath does not affect the equilibration time Eq. (55). Indeed, it is only the shape of the dissipative kernel that affects thermalization. We also note that although the increase of the average energy per oscillator is fairly negligible, since the energy lost by the heavy quark during thermalization is distributed to N oscillators where N is large, thermalization produces entropy because the information stored in the initial conditions becomes irrelevant after the thermalization time [62].

We can combine the results found in this section with those of Section 2.4 to give a rough estimate of the kinematic regime in which the motion of the heavy quark is dominated by diffusion. Indeed, the drag would correspond to a shift of the initial heavy quark momentum by an amount

where $ \gamma = 1/\tau_{\rm eq} $. On the other hand, diffusion amounts to a dispersion of the kinetic energy around its initial value, Eq. (36). Using the fluctuation-dissipation theorem in the form of Eq. (39), which is valid both for the bath in thermal equilibrium and for the CGC-like bath, we can write

Taking the ratio of the last two equations we find

Loosely speaking, the above equation shows that the drag can be neglected with respect to diffusion as long as the initial kinetic energy of the heavy quark is small in comparison with the average energy of the bath. For a bath in thermal equilibrium this average energy is just the temperature. For a CGC-like bath, we have from Eq. (9) $ 2{\cal E} \approx {\cal M} \approx Q_s $ , where $ Q_s $ is the saturation scale.

A purely diffusive process has been observed in [48], where the propagation of heavy quarks in the evolving glasma has been studied neglecting the effect of this propagation on the background gluon field. This is called probe approximation, since the color current carried by the heavy quark is supposed not to affect the gluon field. From the model discussed in this work we can interpret the lack of back-reaction as the lack of a drag coefficient in the equations of motion of the heavy quarks. While we recognize the need for a full numerical solution of the problem of propagation of heavy quarks in the evolving glasma, we can use the present model, and in particular the estimate (58), to state that the pure diffusive approximation should work well as long as the initial kinetic energy of the heavy quark is small in comparison with the saturation scale. For example, taking $ Q_s \approx 2 $ GeV and an initial momentum $ P(0) = 2 $ GeV, which is a fair estimate of the initial average transverse momentum of heavy quarks produced via perturbative QCD [48], we read from Eq. (58) that the ratio between drag and diffusion terms is approximately $ 0.67 $ for the c quark, and $ \approx 0.22 $ for the b quark, meaning that the pure diffusive approximation in this example works pretty well for the b quark but is marginal for the c quark. For higher momenta the back-reaction needs to be considered in order to have a correct description of the evolution of the heavy quark in a dense gluon system.

In this section we analyze the solution of the equation of motion of P at very large time. This is important because it gives information on the equilibrium state of the heavy quark in the bath of oscillators. Moreover, it is useful to understand how the non-Markov dissipative kernel, as well as the out-of-equilibrium oscillators bath, affect the large time state of the heavy quark. From Eqs. (54) and (51) we note that only the poles $ s = \pm i\omega_k $ give a contribution at large times since the other poles lead to an exponential decay of the solution. Therefore, for very large time we can write

and taking into account only the relevant poles we get

We note that although the poles corresponding to $ s+\Gamma(s) = 0 $ do not contribute, the dissipative kernel still appears explicitly in Eq. (60).

We use Eq. (60) to compute the average kinetic energy of the heavy quark. In order to do so we need to square the right hand side of the equation, then take the ensemble average with the distributions (4) or (7), and then sum over the oscillators. This is a pretty straightforward procedure which is not enlightening, therefore we skip the few mathematical steps and write down the final result:

where we have used $ \Gamma = \Gamma(i\omega_k) $, $ \Gamma^* = \Gamma(-i\omega_k) $ and

We briefly note that the second term in the right hand side of Eq. (61) would contribute only if $ \langle\alpha_k^2\rangle \neq \langle\beta_k^2\rangle $. The evaluation of this term turns out to be very simple in the case of a Markov kernel in which $ \Gamma(s) $ is a constant. We will not comment further on this term in the present work, although in principle it would be interesting to study the situation where $ \langle\alpha_k^2\rangle \neq \langle\beta_k^2\rangle $ , since this would correspond to a different kind of an out-of-equilibrium bath in which the different degrees of freedom thermalize at different temperatures.

For the kind of baths considered in this work we have $ \langle\alpha_k^2\rangle = \langle\beta_k^2\rangle = g^2{\cal E}/(m\omega_k^2) $ , with $ {\cal E} = T $ for the system in thermal equilibrium, and $ {\cal E} = {\cal T}_{\rm eff} $ for the CGC-like distribution. Therefore, only the first term on the right hand side of Eq. (61) gives a contribution:

We were able to evaluate the above integral only in the limit of large a, when we take the asymptotic expansion of the integrand up to the order $ 1/a $:

The above equation shows that for a kernel without memory, the heavy quark equilibrates to the average kinetic energy of the bath. The presence of memory, which corresponds to the term $ \propto 1/a $ in Eq. (64), leads to a higher value of the average kinetic energy of the heavy quark if the coupling $ \alpha $ is kept fixed. Clearly, this does not correspond to a measurable increase of the average energy of the bath. The heavy quark loses part of its initial energy, $ \Delta E $, during the thermalization process, and this $ \Delta E $ is distributed to the bath so that the average energy of each oscillator increases by $ \Delta E/N $ , which is negligible in the limit $ N\rightarrow\infty $ . Nevertheless, the entropy increases during thermalization because the information of the initial conditions is lost, as already specified above.

The result in Eq. (64) shows that in presence of the memory the heavy quark does not thermalize to the temperature of the bath: this is expected in our model since memory is always present during evolution, therefore the interaction with the bath is constantly dirtied by the interaction of the heavy quark with its own perturbation and this will be present also at equilibrium. We note that Eq. (64) has been derived assuming that memory is always present: if we make the assumption that our parameter a changes with time we can implement a memory kernel which is Gaussian at the initial time, then becomes a $ \delta $-function at larger times. Doing so results in a perfect equilibration of the heavy quark with the bath.

We have considered a simple model for the diffusion of heavy quarks in a hot bath, the latter being modeled by an ensemble of harmonic oscillators distributed according to either a thermal distribution or to an out-of-equilibrium CGC-like distribution with a saturation scale. The use of the latter is motivated by the recent increased interest of studying the propagation of heavy quarks in a dense medium of out-of-equilibrium gluons that can be produced in high energy nuclear collisions. Therefore, this model might help to shed light on the involved dynamics of the early stages of these collisions avoiding the complications of solving the non-abelian equations of motion of the gluon background coupled to the heavy quarks. Of course, we do not pretend that this model can replace the study of the full problem. We are well aware that a full implementation of the gluon + quarks dynamics is essential to make quantitative predictions. Indeed. this is what we have started in [48]. Nevertheless, the simple model studied here can help to clarify at least the qualitative picture, illustrating the roles in the interplay between the heavy quarks and the bath, and hence suggesting how to improve the present calculations in order to get a more complete picture. The diffusion model is well known in literature since it has been studied in the context of dissipation in quantum mechanics [50–53].

In this model it is possible to study the process of thermalization of heavy quarks, in particular how it depends on the specific form of the distribution used for the bath, as well as on the coupling and the heavy quark mass. Moreover, it is possible to include memory effects in the dissipative kernel, $ \gamma(t-t^{\prime}) $, and to study how it affects thermalization. The memory can be implemented easily by choosing properly the distribution of oscillators, $ {\rm d}N/{\rm d}\omega $, which appears in the very definition of $ \gamma(t-t^{\prime}) $. For the model considered here, a Markov kernel $ \gamma(t-t^{\prime}) = 2\gamma\delta(t-t^{\prime}) $ can be obtained by assuming $ {\rm d}N/{\rm d}\omega\propto\omega^2 $. We studied the memory effects by introducing continuous deformations of the Markov kernel via a Gaussian distribution, which helps to keep the mathematical expressions simple, and at the same time allows for a smooth transition from a kernel without to a kernel with memory. Therefore, the two main merits of this model are the ability to implement memory effects and an out-of-equilibrium dense bath of oscillators.

We remarked that in order to obtain the dissipative term in the equation of motion of the heavy quark, namely the drift, it is necessary to include the back-reaction of the oscillators on the bath. Indeed, from the dynamical equation it is clear that besides the coupling of the quark to the unperturbed oscillators, which is responsible for diffusion, there is also coupling of the heavy quark with the perturbation of the bath induced by the latter, which is responsible for the dissipation. Moreover, we have proved that with the assumption $ \langle p_{0k}^2\rangle = m^2\omega_k^2\langle x_{0k}^2\rangle $ , the fluctuation-dissipation theorem is satisfied both for the bath at thermal equilibrium and for the CGC-like bath, regardless of the shape of the distribution of oscillators.

The solution of the equation of motion for the heavy quark momentum, P, can be obtained easily. Using this solution, we defined the equilibration time as the time necessary for the system to lose information about its initial condition. It was possible to compute analytically the thermalization time by studying the distribution of the zeros of $ s+\Gamma(s) $ in the complex plane, where $ \Gamma(s) $ denotes the Laplace transform of the dissipative kernel. We selected the zero with the smallest negative real part, $ -\chi $, and defined the inverse of the thermalization time as $ 1/\tau_{\rm eq} = \chi $. We found that $ {\rm d}N/{\rm d}\omega $ is the only information about the bath that affects the thermalization time: the distribution of the degrees of freedom $ \langle p_{0k}^2\rangle $ and $ \langle x_{0k}^2\rangle $ does not affect the thermalization time. We also found that keeping the strength of the effective coupling fixed and deforming $ {\rm d}N/{\rm d}\omega $ in order to have a kernel with memory, the thermalization time decreases. We understood this effect as a consequence of the larger density of oscillators, which in turn results in a more efficient way to dissipate energy. We remark that we obtained this result using a specific analytical form of $ {\rm d}N/{\rm d}\omega $, namely a Gaussian. At this moment we are unable to make a general statement on the effect of memory on the thermalization time. It would be interesting to study other kernels with memory to check whether our finding is specific to the Gaussian or is more general, even outside the context of heavy quarks dynamics.

The work presented was partly inspired by [48] , in which we studied the diffusion of heavy quarks in an evolving glasma background. In [48] , we used the probe approximation, which amounts to neglecting the back-reaction of the colored current carried by the heavy quarks on the gluon background, and diffusion-without-drag was found to be the main characteristics of the propagation of heavy quarks in this medium. We were able to understand why in [48] the lack of drift was observed. Indeed, from the equation of motion for the momentum of the heavy quark we noted that dissipation arises from the coupling of the heavy quark to the perturbation induced by the quark itself on the bath of oscillators. Therefore, in order to observe dissipation, the implementation of the back-reaction on the gluon background seems necessary. Nevertheless, we can use this model to estimate roughly the kinematic regime in which the diffusion-without-drag is a fair approximation. Indeed, from the fluctuation-dissipation theorem that we derived for the CGC-like distribution we found that the relative shift of the kinetic energy due to drag and diffusion is just given by $ \approx K_0/{\cal E} $, where $ K_0 $ denotes the initial (non-relativistic) kinetic energy of the heavy quark, and $ {\cal E} $ is the energy-per-oscillator of the bath. For a system in thermal equilibrium $ {\cal E} = T $, while for the CGC-like distribution we can take $ {\cal E} \approx Q_s $, where $ Q_s $ is the saturation scale. Therefore, our model suggests that as long as the initial energy of the heavy quark is smaller than the average energy of the bath, the pure diffusive motion is not a bad approximation.

There are numerous possible studies that can be done in the future. As mentioned above, it would be interesting to explore the different shapes of the dissipative kernel. It also seems interesting to investigate the medium-induced interaction between heavy quarks that might affect the thermalization time, as well as observables like particle-particle correlations. Besides, it is of a certain interest to consider more general baths in which $ \langle p_{0k}^2\rangle \neq m^2\omega_k^2\langle x_{0k}^2\rangle $, which implies that different degrees of freedom of the bath equilibrate at different temperatures. Therefore, an interesting point to study would be to include self-interactions in the bath in order to study its equilibration, besides the propagation of the heavy quark in this dynamical bath. In addition to these problems, another aspect that is worth investigating is the correlation between different degrees of freedom of the bath. As an application of the present work, we intend to explore in a forthcoming article the impact of memory on heavy quark dynamics and its thermalization time at RHIC and LHC energies, which was ignored in earlier calculations.

The authors acknowledge Gabriele Coci, Vincenzo Greco, Lucia Oliva, John Petrucci, Salvatore Plumari and in particular Xingbo Zhao for inspiration, useful discussions and comments on the first version of this article.