Memory effects on energy loss and diffusion of heavy quarks in the quark-gluon plasma

We study the dynamics of heavy quarks in a thermalized quark-gluon plasma with a time-correlated thermal noise, $\eta$. In this case it is said that $\eta$ has memory. We use an integro-differential Langevin equation in which the memory enters via the thermal noise and the dissipative force. We assume that the time correlations of the noise decay exponentially on a time scale, $\tau$, that we treat as a free parameter. We compute the effects of $\tau\neq 0$ on the thermalization time of the heavy quarks, on their momentum broadening and on the nuclear modification factor. We find that overall memory slows down the momentum evolution of heavy quarks: in fact, transverse momentum broadening and the formation of $R_{AA}$ are slowed down by memory and the thermalization time of the heavy quarks become larger. The potential impact on other observables is discussed briefly.


I. INTRODUCTION
A hot and dense phase of nuclear matter, the quarkgluon plasma (QGP), is expected to form in the ultrarelativistic heavy-ion collisions at Relativistic Heavy-ion Collider (RHIC) and the Large Hadron Collider (LHC) energies. Probing and characterizing the bulk properties of QGP is a field of high contemporary interest. Heavy quarks (HQs) [1][2][3][4][5][6][7][8][9][10][11][12] such as charm and beauty are considered as good probes of the system produced in highenergy nuclear collisions. In fact, they are produced in the very early stage due to the hard partonic scatterings on a time scale τ = O(1/m) where m is the rest mass of the quark. Due to their large mass and low abundance, they can propagate in the QGP bringing almost no disturbance to it. Consequently, they act as good probes that can experience the whole evolution of the system created in the collisions, from the very early stage up to the hadronization stage.
The standard approach to study the HQ dynamics in the QGP is following their position and momentum evolution by means of the Langevin equations [13][14][15][16][17][18][19][20][21][22][23][24][25][26] (see also [27]) as well as relativistic kinetic theory [9-12, 21, 28-31]. In the approaches based on the Langevin equation, the thermal noise, η, is usually treated as standard Wiener process, thus without correlations in time. In this work we relax this approximation and analyze the case in which η is time-correlated; this class of stochastic processes is called a process with memory.
The prototype of Langevin equation that we consider in this work is dp dt where p is the momentum of the particle, η is the stochastic term that models the thermal noise, while the integral term on the right hand side is the dissipative force. In previous studies the latter is replaced by −γp where γ is the drag coefficient: this replacement follows from the Fluctuation-Dissipation Theorem (FDT) when η has no time correlations.
Our goal is to analyze the motion of heavy quarks in a quark-gluon plasma, when the correlations of the thermal noise do not decay instantaneously: instead, we assume that these correlations decay over a specific time scale that we call the memory time, τ . Hydrodynamic fluctuations [32,33], diffusion in the evolving Glasma [34][35][36][37][38][39][40][41][42][43], diffusion of electric charge [44], dilepton yields [45] and the electric conductivity of the quark-gluon plasma [46] are some of the physical problems where memory can play a role; for these τ lies in between 0.1 and 10 fm/c. In this study we treat τ as a free parameter and study its effect on a few physical quantities, namely the momentum broadening of heavy quarks and on the nuclear modification factor, R AA . For the sake of simplicity we consider the interaction of heavy quarks with a thermalized quark-gluon plasma at a fixed, constant temperature T ; the diffusion coefficients that we use in the calculations are that obtained by perturbative QCD (pQCD) for high T , and by a quasiparticle model (QPM) for low T , while the dissipative kernel is related to the thermal noise by the Fluctuation-Dissipation Theorem (FDT).
We anticipate the main result, namely that memory delays the dynamics of the heavy quarks in the QGP: we show this by studying the thermalization time, the momentum broadening and the time evolution of R AA . The latter in particular can be potentially of interest for the phenomenology of heavy quarks in the QGP, due to the fact that the slower evolution of R AA would require the use of larger diffusion coefficients in order to reproduce the experimental data and this in turn would require stronger interactions of the heavy quarks with the bulk, arXiv:2203.06712v1 [hep-ph] 13 Mar 2022 potentially leading to a larger v 2 .
The plan of the article is the following: in section II we present the calculations for the equilibration time and the momentum broadening for the Brownian motion with memory in the nonrelativistic limit; in section III we discuss the numerical implementation of the Langevin equation with an integral kernel, while in section IV we present our results. Finally, in section V we draw our conclusions.

II. NONRELATIVISTIC LIMIT
For the sake of illustration, we consider here a simple one-dimensional motion of a heavy particle with mass m in the nonrelativistic limit. Most of the calculations presented here have been obtained firstly in [47], where a gaussian correlator has been considered. Here we consider an exponential correlator instead, therefore we skip many details that have been given in [47], and limit ourselves to write explicit results that stand for the exponential correlator.
The Langevin equation for momentum p is where η is the stochastic term that models the noise, while the integral term on the right hand side is the dissipative force. The formal solution of Eq. (2) can be obtained by means of Laplace transforms, namely where p 0 = p at t = 0, and Γ and Ξ denote the Laplace transforms of the dissipative kernel and of the noise respectively; the integral is understood on a Bromwich contour that leaves all the singularities of the integrand on its left side. We assume that η is a gaussian random variable with correlators given by Moreover, we assume that η represents the thermal noise in a thermalized bath at temperature T . From the FDT we have then In order to analyze the thermalization time and the momentum broadening of the heavy particle we need to evaluate the following averages: where we have put that describes momentum randomization due to the propagation of the particle in the bath, and where L −1 [h](t 1 , t 2 ) is the two-dimensional inverse Laplace transform of h(s 1 , s 2 ) that depends on (t 1 , t 2 ).
In the above equation we have put (11) In this work we consider an exponential correlator in Eqs. (4) and (5), namely which allows to solve the problem analytically. In the following we call τ the memory time, since it sets the time scale over which time correlations of the noise decay. Note that lim τ →0 f (t) = δ(t) in the distributional sense, therefore it is possible to interpolate between the local and the nonlocal kernel by changing the value of τ . The Laplace transform of the dissipation kernel is .
A. Thermalization Firstly, we examine the thermalization of the heavy particle, which consists in the loss of information about the initial condition and in the equilibration of its kinetic energy with the bath, E = T /2 as required by the Equipartition Theorem. From Eqs. (7), (9) and (13), considering that the zeroes of s+Γ(s) are the solutions of the equation s + τ s 2 + γ = 0, we get by a straightforward application of the residues theorem with A ≡ √ 1 − 4γτ . Before moving to the thermalization time we pause on the early, pre-thermalization behavior of p(t) , in order to emphasize the differences between the motion with and without memory. In order to simplify the discussion we assume that τ 1/γ where 1/γ represents the thermalization time for processes without memory. For t τ we get while for τ t 1/γ we get The memory changes the pre-thermalization evolution of p(t) from linear to quadratic; in particular, this implies that the thermalization of the particle in a bath with memory is slower than the one in a bath without memory with the same drag coefficient, γ. This result is confirmed by the calculation of the thermalization time that makes use of the full result (15). We define the thermalization time, τ therm , such that p(τ therm ) = p 0 /e with p(t) given by Eq. (15). The results of this calculation are shown in Fig. 1 where we plot τ therm as a function of the memory time, both measured in units of 1/γ. Thermalization time increases with τ in agreement with the discussion above; the quantitative effect of memory is negligible for γτ 1, but becomes substantial ≈ 20% already for γτ ≈ 0.5.
The qualitative behavior of the thermalization time can be understood in simple terms. In fact, thermalization implies the loss of the correlation with the initial condition: the heavy particle equilibrates with the bath regardless of its initial momentum distribution. If the thermal noise of the bath has memory, the momentum evolution from time t to t + ∆t does not delete totally the information of p(t) because of the correlations of the noise. Therefore, it is natural to expect that the loss of the information about the initalization requires more time.

B. Momentum broadening
Next we turn to the momentum broadening, σ p in Eq. (8). The starting point is the computation of J in Eq. (10) with correlator given by Eq. (14). We get where A has been defined right after Eq. (15). Although Eq. (18) is exact, it is quite cumbersome, therefore it is convenient to analyze a few limit situations in which the result (18) is manageable; after that, we will present the full result (18) computed numerically. As in the previous subsection, for the sake of illustration we assume that τ 1/γ: this assumption will be removed in the full numerical calculation, see Fig. 2. From the results presented in the previous subsection in this regime we have τ therm ≈ 1/γ, see Fig. 1.
The dissipative force becomes relevant on the time scale t ≈ τ therm ≈ 1/γ, thus for tγ 1 we can ignore it and put γ = 0 in Eq. (18). We get In the early stage t τ the above result gives In particular, the result (21) agrees with that we would find for the Brownian motion without memory. We notice that the effect of a finite memory time is to slow down the momentum broadening of the heavy particle, changing the early time evolution from linear to quadratic. For late times t 1/γ, t τ , the memory with the exponential kernel has no effect on the equilibration value of momentum broadening: in fact, in the asymptotic limit γt 1 we get from Eq. (18) in agreement with the standard result for the Brownian motion.

III. LANGEVIN EQUATION WITH MEMORY: NUMERICAL IMPLEMENTATION
In this work, we solve the Langevin equation for the heavy quarks in a bath with a colored noise; the latter has correlations at different times. In order to generate this noise we introduce an ancillary stocastic process, h(t), that evolves simultaneously to (and independently of) the heavy quarks, built up in such a way its correlator at different times does not vanish. In this section we firstly define the ancillary process and specify its correlations at different times; then, we formulate the Langevin equation where the heavy quark is coupled to h(t) and discuss the numerical discretization scheme adopted in the calculations.

A. Ancillary stochastic process
Let us consider the stochastic process a that satisfies the Langevin equation where ξ is a gaussian white noise, a is assumed to be dimensionless, so the 1/α in (25) is put to balance the time dimension carried by the δ−function because ξ is dimensionless too. The formal solution of Eq. (23) is given by where a 0 = a(t = 0). Clearly we have We define the fluctuating field This satisfies Eq. (23) with h(t = 0) = 0, that we rewrite for the sake of future reference: We baptize h as the ancillary process, because we use it as an additional stocastic process to generate the colored noise for the Langevin equation of the heavy quarks, see below.
From the very definition of h it is easy to see that h(t) = 0. Instead the correlator of h at different times is not a δ−function: it can be obtained easily from Eq. (26), namely Therefore, h is a process with memory and it can be used in any Langevin equation. From Eq. (30) it is obvious that asymptotically namely correlations are washed out on the time scale 1/α ≡ τ . Note that for α → +∞ we get, from Eq. (30), that is the process h becomes a standard white noise in the limit τ → 0, as expected. Before going ahead, it is useful to comment on Eq. (30): we note that the correlator is not a function of t−t but of t and t separately. It is convenient to fix t and analyze the correlator for t ≥ t . Numerically Eq. (29) can be discretized in the usual manner by the replacements where ∆t corresponds to the discrete time step implemented in the numerical calculation. With these we have ζ(t) will be implemented as a white noise with variance equal to one.

B. Application to the Langevin equation
Next we turn to the solution of the Langevin equation, (2), for heavy quarks in the relativistic limit. We assume the FDT in the relativistic limit, namely with E = p 2 + M 2 . Moreover, the process η(t) in (2) is assumed to satisfy η = 0 and where τ is the memory time and g is a dimensionless function that defines the correlation of the noise; for simplicity we assume it satisfies g(0) = 1. In the case of a Markov process τ → 0 and g has to satisfy the condition In this work we generate the noise η by means of the ancillary process h introduced in the previous subsection with α = 1/τ . More precisely we assume that see Eq. (30). Therefore, in Eq. (2) we put see also Eq. (37). Using Eqs. (32) and (37) we note that in the limit τ → 0 we get namely we recover the time correlation of the standard Brownian motion. We have noted that the ancillary process (29) develops substantial correlations after an initial transient stage that lasts for t ≈ τ ; after this transient, the correlator is approximately invariant under time translation and decays exponentially on the time scale τ . We call this stage as the exponential decay regime. In the numerical calculations we start the ancillary process and let it run up to some time t 0 , leaving heavy quarks frozen in momentum and coordinate space; then, when the noise is in the exponential decay regime, we start the evolution of the heavy quarks including their interaction with the noise itself. In this regime the correlator takes the form With the rescaling (40) the Langevin equation (2) that we solve for t > t 0 . In (43) we put in agreement with the relativistic form of the FDT. The time-discretized version of Eq. (43) is given by Note that differently from the Markov process, the stochastic term in (45) does not come with the √ ∆t: this is so because now the noise is not a Wiener process due to the memory. We implement a simple iterative scheme to solve (45), namely with s 0 = t 0 , s N −1 = t − ∆t and s N = t. This scheme leads to the self-consistent solution Equations (35) and (45) allow us to implement the momentum evolution of the heavy quark in a bath with a colored noise. h in Eq. (45) is given by the solution of the ancillary Langevin equation (35). This means that at each time step one has to solve (35) [with the initial condition h(t = 0) = 0] and (45) simultaneously. This procedure is different from the one adopted in the literature in the case of a Markov process, in which the white noise would be extracted randomly from a Gaussian distribution at each time step.
For the numerical implementation of the 3-dimensional Langevin equation, we are solving equations (35) and (45) simultaneously for the 3-components (h x ,h y , h z , p x , p y , p z ) to study HQ momentum evolution coupled with the coordinate evolution where dr i is the shift of the coordinate in each time step dt. E and p are the energy and momentum of the heavy quark respectively. The HQ transport coefficients are computed as follows. At high T we model the thermalized bath by a QGP made of massless quarks and gluons, and use the pQCD kinetic coefficients for the processes c → c , where denotes either a massless gluon or a quark in the bath. The diffusion coefficients in this case are well known and can be found in the literature, see for example [49,50]. The squared invariant scattering amplitudes are the Combridge ones, that includes s, u and t channel and their interferences terms. The infrared divergence associated with the t-channel diagrams is screened by the Debye mass, m D = g(T )T . For more details see earlier works [28,51]. On the other hand, at low T we model the bath with a gas of quarks and gluons quasiparticles, using the so-called quasi-particle models (QPMs); in these models the quasiparticle masses are tuned in order to reproduce Lattice QCD thermodynamics [52,53].The QPMs account effectively for the non-perturbative effects for σp versus time for a one-dimensional motion of charm quarks in the nonrelativistic limit, see Eq. (6). σp is defined in Eq. (8). We have considered two values of the memory time, τ . Analytical results correspond to Eq. (18).
T close to the quark-hadron transition temperature, T c . The main feature of the QPM is that the effective coupling is significantly stronger than the one of pQCD near T c , which enhances the HQ-bulk inteactions. We have evaluated the diffusion coefficient within the QPM starting from the effective coupling with massive quarks and gluons. For details we refer to earlier works [28,51].

A. Non-relativistic check
In order to check the discretization scheme of the Langevin equation with correlated noise (45), we plot σ p versus time in Fig. 2; the definition of σ p is given in Eq. (8), and the Langevin equation has been solved for a one-dimensional motion of charm quarks. We have solved the equation in the nonrelativistic limit as explained in section II, that corresponds to replace the kinetic energy with the mass of the quark in the Fluctuation-Dissipation theorem, see Eq. (6); moreover, we have used a constant diffusion coefficient D = 0.2 GeV 2 /fm for illustrative purposes only. Here T referes to the bath temeprature. The green and orange solid lines correspond to τ = 0.2 and τ = 0.01 fm/c. For comparison, in the same figure we show by dashed lines the analytical result (18), that is valid in the nonrelativistic limit: the agreement between the two results is excellent, showing that our numerical scheme works properly. We notice that the memory slows down the evolution of σ p as anticipated in section II. Moreover, for τ = 0.2 fm/c we notice that the diffusive motion is characterized by the initial nonlinear increase of σ p that turns into a linear regime before the drag force becomes relevant and leads to the thermalization of the heavy quark. In the case of the smaller τ the charm quark enters the linear regime immediately, then equilibrates with the medium.

B. Transverse momentum broadening
In this subsection we analyze momentum broadening of heavy quarks in a hot medium; analogously to the previous section we define where p T = p 2 x + p 2 y is the transverse momentum. In Fig. 3 we plot σ p versus time for three temperatures and two values of τ , namely τ = 0.2 fm/c (solid lines) and τ = 1 fm/c (broken lines). Calculations correspond to the cham quark. For the three cases considered it is clear that increasing the memory time results in the slowing down of momentum broadening; the effect is more visible at large temperature, where the drag force is less effective. At small temperature it is still possible to measure some difference between the results with the two memory times in the early evolution, then for time t ≥ 5 fm/c the results with and without memory coincide.
In Fig. 4 we plot a selection of the results shown in Fig. 3 for τ = 0.2 fm/c, zooming on the early time evo- lution of σ p . We notice the nonlinear increase of σ p , induced by memory in agreement with the discussion of section II, followed by a linear enhancement. In this regime the charm quarks experience an almost diffusive motion, in the sense that the energy loss due to the drag force is still negligible. Qualitatively the different regimes appear also for the small temperature case in the figure: however, in this case the drag force is stronger so the linear regime lasts for a shorter fraction of the evolution, then σ p bends and eventually saturates, signaling the equilibration of the charm with the medium.

C. Nuclear modification factor
In this subection we analyze the modification factor R AA , defined as where (dN/d 2 p T ) t denotes the spectrum of charm quarks at time t and (dN/d 2 p T ) FONLL denotes the spectrum at the initialization time, To this end, at the formation time we assume the prompt spectrum obtained within Fixed Order + Next-to-Leading Log (FONLL) QCD that reproduces the D-mesons spectra in pp collisions after fragmentation [54,55]. for charm quarks; the slope of the spectrum has been calibrated to a collision at √ s = 5.02 TeV. R AA (p T ) = 1 implies that charm quarks experience interactions with the gluon medium, causing a change in their spectrum. Our motivation is to highlights the impact of memory on R AA (p T ).
In Fig. 5 we plot R AA versus p T , for two values of the memory time τ . The shape of R AA for T = 1 GeV and up to p T ≈ 5 GeV is in agreement with the one advertised in [34,35,40]: this is the result of the diffusion-dominated propagation of the heavy quarks in the hot medium at high temeprature, that effectively diffuses low momentum charm quarks to higher momentum states. This explains why R AA remains smaller than one up to p T ≈ 2 GeV. For larger p T the energy loss is important and particles migrate to lower momentum states. At low temperature the drag force is dominant in the whole range of p T , therefore the general tendency is that particles move to lower p T states so R AA stays greater than one up to p T ≈ 2 GeV.
The effect of τ = 0 is to slow down the formation of R AA . We notice a sizable impact of memory on the R AA . To make this poing clearer, in Fig. 6 we show R AA for τ = 0 and τ = 1 fm/c at T = 0.25 GeV obtained within QPM. We notice the slower evolution of R AA when τ = 0.
In Fig. 7 we plot R AA at t = 3 fm/c, for τ = 0 and τ = 1 fm/c, at T = 0.25 GeV. Solid lines correspond to the results of the QPM model, while dashed lines stand for the pQCD calculations. We notice that R AA obtained within the QPM model differs considerably from the results of pQCD, which is in agreement with expectations because the cross sections within the former are enhanced with respect to those of the latter, implying that diffusion and drag coefficients are larger.
The results in Fig. 7 show that memory slows down the dynamics of the heavy quarks in the whole p T range. In fact, R AA at large p T in the process with memory stays above that without memory, meaning that the large p T particles have lost less energy during the evolution. In other words, memory slows down the energy loss. Ergo, also thermalization in presence of a time-correlated noise is retarded. Similarly, the diffusion to low p T is slower when memory is present, since R AA in this case remains lower than the one without memory. These effects are more evident in the case of the QPM while they are smaller (but still present) in the pQCD case.
The results discussed in this section allow us to discuss a potential impact that memory might have on more sophisticated phenomenological calculations. As a matter of fact, our results can be summarized by stating that memory slows down the heavy quark dynamics. For example, the formation of R AA of heavy quarks is slower if the thermal fluctuations of the bath are time correlated. Therefore, if the diffusion coefficient is tuned in order to reproduce the experimental R AA , then a larger coefficient is needed when the bath has memory. The larger diffusion coefficient will then enhance the elliptic flow, v 2 . The memory effect has the potential to alter the heavy quark R AA -v 2 dynamics and can improve the simulations description of heavy quark R AA -v 2 [51]. Memory effect may affect other observables like the heavy quark directed flow, particle correlations and so on. Within our present study it is impossible to predict quantitatively this change due to our simplifying assumptions; these simplifications are partly justified by the fact that it is the first time, to our knowledge, that the effects of time correlated fluctuations on heavy quark observables are computed. A more detailed, quantitative study will be the subject of future investigations.

V. CONCLUSIONS
We have studied the effects of a time-correlated thermal noise of a thermalized quark-gluon plasma on the energy loss and the diffusion of heavy quarks. In this case the time correlation of the thermal noise does not decay instantaneously, as instead it is assumed for the standard Brownian motion. We have considered a simple situation in which the time correlations of the noise decay exponentially on the time scale τ , called the memory time, and treated τ as a free parameter. We have considered an integro-differential Langevin equation for the heavy quark momentum, taking into account memory effects both in the thermal noise and in the dissipative force.
We have computed several quantities that characterize the dynamics of heavy quarks in the bath, namely the thermalization time and the transverse momentum broadening. Then we turned to the nuclear modification factor, R AA , that we have computed using kinetic coefficients from pQCD and quasiparticle models. Our results suggest that memory delays the dynamics of the heavy quarks in the QGP: indeed, memory slows down momentum broadening as well as the formation of R AA , retards the energy loss and thus increases the thermalization time.
Our work could be of interest for the phenomenology of heavy quarks in the QGP: in fact, the slower evolution of R AA would require the use of larger diffusion coefficients in phenomenological calculations, for reproducing the experimental R AA and this in turn would require stronger interactions of the heavy quarks with the bulk, potentially leading to a larger v 2 . The investigation of this point requires a refined study with realistic initial conditions, including the initial geometry of the fireball, and will be matter for future studies. In addition to this problem, it is of a certain interest to analyze processes with a long tail memory, in which the time correlations of the noise do not decay exponentially but as power laws. These processes can be generated via Langevin equations with fractional derivatives [56,57]; this problem will also be the subject of future studies.