-
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$ N = \int_{-\infty}^\infty\frac{{\rm d}N}{{\rm d}\omega}{\rm d}\omega. $
(1) 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]$ H = \frac{P^2}{2M} +V(Q)+ \sum\limits_{k = 1}^N \left( \frac{p_k^2}{2m} + \frac{m\omega_k^2 x_k^2}{2} \right) -g Q\sum\limits_{k = 1}^N x_k + {\cal H}_R, $
(2) 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$ \sum\limits_{k = 1}^N f(\omega_k)\rightarrow \int_{-\infty}^{+\infty}\frac{{\rm d}N}{{\rm d}\omega}f(\omega) {\rm d}\omega. $
(3) For a bath in thermal equilibrium, the initial conditions for the oscillators are distributed randomly according to the standard Boltzmann distribution at temperature T,
$ f(p_{0k},x_{0k}) = {\rm e}^{-E^p_k/T}{\rm e}^{-E^x_k/T}; $
(4) where we have used
$ E^p_k = \frac{p_{0k}^2}{2m}, $
(5) $ E^x_k = \frac{m\omega^2_k x_{0k}^2}{2}, $
(6) 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$ \frac{1}{1+{\rm e}^{(E^p_k - {\cal M})/{\cal T}}}\frac{1}{1+{\rm e}^{(E^x_k - {\cal M})/{\cal T}}}, $
(7) 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,$ {\cal T}_{\rm eff} = {\cal T}\frac{{\rm Li}_{3/2}(-z)}{{\rm Li}_{1/2}(-z)}, $
(8) 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$ {\cal T}_{\rm eff} = \frac{2{\cal M}}{3}. $
(9) 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
$ \frac{{\rm d}^2 x_k}{{\rm d}t} + \omega_k^2 x_k = \frac{g}{ m}Q. $
(10) This equation can be easily solved exactly using the retarded Green's function, namely
$ x_k(t) = x_k^h(t) - \frac{g}{ m} \int_0^t G_R^k(t-t^{\prime}) Q(t^{\prime}) {\rm d} t^{\prime}, $
(11) where we have defined
$ {\cal G}_R^k(x) = \theta(x)G_R^k(x), $
(12) with
$ {\cal G}_R^k(x) $ denoting the retarded Green's function of the oscillator,$ \left(\frac{{\rm d}^2}{{\rm d}t^2} + \omega_k^2\right){\cal G}_R^k(t-t^{\prime}) = -\delta(t-t^{\prime}).$
(13) Moreover, we have defined
$ x_k^h(t) = \alpha_{k}\cos\omega_k t + \beta_k \sin\omega_k t, $
(14) 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
$ \frac{{\rm d}P}{{\rm d}t} = -\frac{\partial V}{\partial Q} -\frac{\partial{\cal H}_R}{\partial Q}+ g\sum\limits_k x_k(t). $
(15) Defining the source term,
$ \xi(t) $ , as$ \xi(t) = g\sum\limits_{k = 1}^N x_k^h(t), $
(16) 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$ \frac{{\rm d}P}{{\rm d}t} +\frac{ g^2}{ m}\sum\limits_{k = 1}^N \int_0^t G_R(t-t^{\prime}) Q(t^{\prime}) {\rm d} t^{\prime} = -\frac{\partial V}{\partial Q} -\frac{\partial{\cal H}_R}{\partial Q} + \xi(t). $
(17) 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
$ \begin{split} \int_0^t G_R(t-t^{\prime}) Q(t^{\prime}) {\rm d} t^{\prime} =& -F(0) Q(t) + F(t) Q(0) \\&+ \frac{1}{M}\int_0^t F(t-t^{\prime}) P(t^{\prime}) {\rm d} t^{\prime}, \end{split}$
(18) 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$ \begin{split} \frac{{\rm d}P}{{\rm d}t} + &\frac{g^2}{mM}\sum\limits_{k = 1}^N \int_0^t F(t-t^{\prime}) P(t^{\prime}) {\rm d} t^{\prime} = -\frac{\partial V}{\partial Q} -\frac{\partial{\cal H}_R}{\partial Q} \\&+ \sum\limits_{k = 1}^N F(0) Q(t) - \sum\limits_{k = 1}^N F(t) Q(0)+ \xi(t). \end{split}$
(19) 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$ {\cal H}_R = \frac{1}{2}\sum\limits_k F(0) Q^2 $
(20) 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$ \frac{{\rm d}P}{{\rm d}t} + \frac{g^2}{mM}\sum\limits_{k = 1}^N \int_0^t F(t-t^{\prime}) P(t^{\prime}) {\rm d} t^{\prime} = -\frac{\partial V}{\partial Q} - F(t) Q(0)+ \xi(t). $
(21) 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$ \frac{{\rm d}P}{{\rm d}t} + \frac{g^2m}{M}\sum\limits_{k = 1}^N \int_0^t F(t-t^{\prime}) P(t^{\prime}) {\rm d} t^{\prime} = \xi(t). $
(22) 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,$ G_R^k(t-t^{\prime}) = -\frac{\sin[\omega_k(t-t^{\prime})]}{\omega_k}; $
(23) the primitive of
$ G_R^k $ is easily calculated. Defining the dissipative kernel as$ \gamma(t-t^{\prime}) = \frac{g^2 }{mM} \sum\limits_{k = 1}^N \frac{\cos[\omega_k(t-t^{\prime})]}{\omega_k^2}. $
(24) we can write the dynamical equation for P as
$ \frac{{\rm d}P}{{\rm d}t} + \int_0^t \gamma(t-t^{\prime}) P(t^{\prime}){\rm d}t^{\prime} = \xi(t). $
(25) 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
$ \gamma(t-t^{\prime}) = \frac{g^2 }{2mM} \int_{-\infty}^\infty {\rm d} \omega \left(\frac{1}{\omega^2}\frac{{\rm d}N}{{\rm d}\omega}\right) {\rm e}^{{\rm i}\omega(t-t\prime)}, $
(26) 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$ \frac{{\rm d}P}{{\rm d}t} + \gamma P = \xi(t). $
(27) 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$ \begin{split} \langle\xi(t)\xi(t^\prime)\rangle =& \frac{1}{2}\sum\limits_k\left\{ \langle \alpha_{k}^2\rangle + \langle \beta_{k}^2\rangle\right\}\cos[\omega_k(t-t^\prime)] \\ &+\frac{1}{2}\sum\limits_k\left\{ \langle \alpha_{k}^2\rangle - \langle \beta_{k}^2\rangle\right\}\cos[\omega_k(t+t^\prime)]. \end{split} $
(28) The above equation is general. If we now make the further assumption that
$ \langle p_{0k}^2\rangle = m^2\omega_k^2\langle x_{0k}^2\rangle, $
(29) we get
$ \langle\xi(t)\xi(t^{\prime})\rangle = \frac{g^2}{m^2}\sum\limits_k \frac{\langle p_{0k}^2\rangle}{\omega_k^2} \cos[\omega_k(t-t^{\prime})]. $
(30) Assuming that
$ \langle p_{0k}^2\rangle $ is independent of the oscillators, which happens for the distributions used in this work, we can write$ \langle\xi(t)\xi(t^{\prime})\rangle = \frac{g^2\langle p_{0}^2\rangle}{m^2}\sum\limits_{k = 1}^N \frac{\cos[\omega_k(t-t^{\prime})]}{\omega_k^2} , $
(31) 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
$ \gamma(t-t^{\prime}) = \frac{m}{M\langle p_{0}^2\rangle}\langle\xi(t)\xi(t^{\prime})\rangle. $
(32) 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
$ P(t) = P(0) + \int_0^t {\rm d}t^{\prime}\xi(t^{\prime}).$
(33) Squaring and taking the ensemble average leads to
$ \langle P(t)\rangle = P(0), $
(34) $ \left\langle P(t)^2 \right\rangle = P(0)^2 + \frac{2 g^2 {\cal E}}{m} \int {\rm d}\omega\frac{{\rm d}N}{{\rm d}\omega}\frac{1-\cos\omega t}{\omega^4}, $
(35) 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$ \left\langle \frac{P(t)^2}{2M} \right\rangle - \frac{P(0)^2}{2M} = 2Dt, $
(36) where we have defined the diffusion coefficient of kinetic energy as
$ D = \frac{ \pi g^2 {\cal E}}{2mM}\left(\frac{1}{\omega^2}\frac{{\rm d}N}{{\rm d}\omega}\right). $
(37) 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
$ \gamma = \frac{\pi g^2}{2mM}\left(\frac{1}{\omega^2}\frac{{\rm d}N}{{\rm d}\omega}\right), $
(38) and comparing this with Eq. (37), we find
$ \frac{D}{\gamma} = {\cal E}; $
(39) 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$ \frac{{\rm d}N}{{\rm d}\omega} = \frac{2N}{a^3\sqrt{\pi}}\omega^2 {\rm e}^{-\omega^2/a^2}, $
(40) 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)
$ \gamma(t-t^{\prime}) = \frac{2\alpha^2}{M}\Phi_a(t-t^{\prime}), $
(41) where
$ \Phi_a(t-t^{\prime}) = \frac{a}{2\sqrt{\pi}}\exp\left(-\frac{a^2(t-t^{\prime})^2}{4}\right), $
(42) and
$ \alpha^2 = \frac{g^2 N }{a^3 m}\sqrt{\pi}. $
(43) 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
$ {\cal L}_s(f) = F(s) = \int_0^\infty {\rm e}^{-st}f(t){\rm d}t.$
(44) It is then straightforward to write Eq. (25) as
$ s{\tilde P}(s)-P(0)+\Gamma(s){\tilde P}(s) = \Xi(s). $
(45) 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$ P(t) = \frac{1}{2\pi i}\int_{\sigma-i\infty}^{\sigma+i\infty}\frac{\Xi(s)+P(0)}{s+\Gamma(s)}{\rm e}^{st}{\rm d}s, $
(46) 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:
$ \begin{split} P(t) =& \frac{P(0)}{2\pi i}\int_{\sigma-i\infty}^{\sigma+i\infty}\frac{1}{s+\Gamma(s)}{\rm e}^{st}{\rm d}s \\&+ \frac{1}{2\pi i}\int_{\sigma-i\infty}^{\sigma+i\infty}\frac{\Xi(s)}{s+\Gamma(s)}{\rm e}^{st}{\rm d}s . \end{split} $
(47) 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
$ \Gamma(s) = \frac{\alpha^2}{M} {\rm e}^{\frac{s^2}{a^2}}E_a\left(s/a\right), $
(48) where N is the number of oscillators in the bath and we have used
$ E_a(x) = \frac{2}{\sqrt{\pi}}\int_x^\infty {\rm e}^{-t^2}{\rm d}t.$
(49) 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$ \Xi(s) = \sum\limits_{k = 1}^N\left(\frac{\alpha_ks+\beta_k\omega_k}{s^2+\omega_k^2}\right). $
(50) Therefore, we can write the exact solution in the form
$ \begin{split} P(t) =& \frac{P(0)}{2\pi i}\int_{\sigma-i\infty}^{\sigma+i\infty}\frac{1}{s+\Gamma(s)}{\rm e}^{st}{\rm d}s \\&+ \frac{1}{2\pi i}\sum\limits_{k = 1}^N\int_{\sigma-i\infty}^{\sigma+i\infty} \left(\frac{\alpha_k s+\beta_k\omega_k}{s^2+\omega_k^2}\right)\frac{1}{s+\Gamma(s)}{\rm e}^{st}{\rm d}s. \end{split}$
(51) 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$ \chi^1 = \frac{\alpha^2 }{M} +\frac{2}{a\sqrt{\pi}}\frac{\alpha^4}{M^2} + O\left(\frac{\alpha^6}{M^3}\right), $
(52) $ \kappa_r^n = na + \frac{\alpha^2 }{M} +\frac{2}{a\sqrt{\pi}}\frac{\alpha^4}{M^2} + O\left(\frac{\alpha^6}{M^3}\right),\; \; \; n = 1,2,\dots $
(53) 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$ \begin{split} P(t) = & P(0) {\rm e}^{-\chi^1 t} + {\rm e}^{-\chi^1 t} \sum\limits_{k = 1}^N \frac{-\alpha_k\chi^1+\beta_k\omega_k}{(\chi^1)^2+\omega_k^2} \\ &+\sum\limits_{k = 1}^N \left[\frac{i\alpha_k\omega_k+\beta_k\omega_k}{2i\omega_k}\frac{{\rm e}^{{\rm i}\omega_kt}}{i\omega_k+\Gamma(i\omega_k)} + {\rm c.c.}\right], \end{split} $
(54) 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$ \tau_{\rm eq} = \frac{M}{\alpha^2} \left(1 - \frac{2\alpha^2}{a \sqrt{\pi} M } + O(\alpha^4/a^2 )\right). $
(55) 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
$ |P(t)^2 -P(0)^2|_{\rm drag} \approx 2P(0)^2 \gamma t, $
(56) 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$ |P(t)^2 -P(0)^2|_{\rm diffusion} = 4M\gamma{\cal E}t. $
(57) Taking the ratio of the last two equations we find
$ \frac{|P(t)^2 -P(0)^2|_{\rm drag}}{|P(t)^2 -P(0)^2|_{\rm diffusion}} \approx \frac{1}{{\cal E}}\frac{P(0)^2}{2M}. $
(58) 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$ P(t) = \frac{1}{2\pi i}\sum\limits_{k = 1}^N\int_{\sigma-i\infty}^{\sigma+i\infty} \left(\frac{\alpha_k s+\beta_k\omega_k}{s^2+\omega_k^2}\right)\frac{1}{s+\Gamma(s)}{\rm e}^{st}{\rm d}s, $
(59) and taking into account only the relevant poles we get
$ P(t) = \sum\limits_{k = 1}^N \left(\frac{i\alpha_k +\beta_k }{2i}\frac{{\rm e}^{{\rm i}\omega_kt}}{i\omega_k+\Gamma(i\omega_k)} + {\rm c.c.}\right). $
(60) 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:
$ \begin{split} \langle P^2\rangle_\infty =& \frac{1}{2}\sum\limits_k \frac{\langle\alpha_k^2\rangle + \langle\beta_k^2\rangle} {(\Gamma^\star - i\omega_k)(\Gamma + i\omega_k)}\\ &+ \frac{1}{4}\sum\limits_k {\cal F}_k(t) \frac{(\langle\alpha_k^2\rangle - \langle\beta_k^2\rangle) } {(\Gamma^\star - i\omega_k)^2(\Gamma + i\omega_k)^2}, \end{split}$
(61) where we have used
$ \Gamma = \Gamma(i\omega_k) $ ,$ \Gamma^* = \Gamma(-i\omega_k) $ and$ {\cal F}_k(t) = {\rm e}^{-2{\rm i}\omega_k t}(\Gamma + i\omega_k)^2 + {\rm e}^{+2{\rm i}\omega_k t}(\Gamma^\star - i\omega_k)^2. $
(62) 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:$ \langle P^2\rangle_\infty = \frac{g^2{\cal E}}{ m} \int_{-\infty}^{+\infty} \frac{1}{\omega^2}\frac{{\rm d}N}{{\rm d}\omega} \frac{{\rm d}\omega} {(\Gamma^\star - i\omega)(\Gamma + i\omega)}. $
(63) 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 $ :$ \left\langle \frac{P^2}{2M}\right\rangle_\infty = \frac{{\cal E}}{2} \left[ 1+\frac{\alpha^2}{aM\sqrt{\pi}} \right]. $
(64) 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.
Classical model for diffusion and thermalization of heavy quarks in a hot medium: memory and out-of-equilibrium effects
- Received Date: 2019-05-23
- Available Online: 2019-09-01
Abstract: We consider a simple model for the diffusion of heavy quarks in a hot bath, modeling the latter by an ensemble of oscillators distributed according to either a thermal distribution or to an out-of-equilibrium distribution with a saturation scale. In this model it is easy to introduce memory effects by changing the distribution of oscillators: we model them by introducing a Gaussian distribution,