Evolution of heavy quark distribution function in quark-gluon plasma: using the Iterative Laplace Transform Method

The"Iterative Laplace Transform Method"is used to solve the Fokker-Planck equation for finding the time evolution of the heavy quarks distribution functions such as charm and bottom in quark gluon plasma. These solutions will lead us to calculation of nuclear suppression factor RAA. The results have good agreement with available experiment data from the PHENIX collaboration.


I. Introduction
Colliding heavy ions in Relativistic HeavyIion Collider (RHIC) and Large hadron Collider(LHC) experiments generate temperatures many thousand times the temperature of the sun. These collisions are recreating in the laboratory conditions similar to those just after the big bang. Under these extreme conditions, protons and neutrons melt, freeing the quarks from their bonds with the gluons. This new phase of matter is called quark-gluon plasma. This highly excited state of matter displays properties similar to a nearly perfect fluid, which has been successfully described by hydrodynamic models [1,2,3,4].This medium is probed by hard processes with scales ranging from a few GeV to several of T eV . The existence of such a phase and its properties are key issues in the theory of quantum chromodynamics (QCD). Basically, heavy partons, which are produced in early stage of heavy ion collisions, have relatively high density. Such heavy partons could induce a large amount of energy loss for hard patrons produced in the initial stage of the collision. So, the observed jet quenching in such collisions can indicate the formation of this strongly interacting dense matter [5,6,7,8,9,10,11]. Charm and bottom quarks play crucial role in this endeavor because they are the witness to the entire space-time evolution of the system and due to their large masses, a memory of their interaction history may be preserved [12,13,14,15].
Their mass is significantly larger than the temperature of the plasma they are formed in, M ≫ T c . They are produced in the early stage of heavy ion collisions and they do not dictate the bulk properties of the matter. In addition the heavy quark thermalization time in perturbative QCD is estimated to be of order of 10 − 15 fm/c and 25 − 30 fm/c

II.The heavy quark momentum evolution in FPE dynamics
The Quark Gluon Plasma (QGP) is formed at the time τ i after the impact. Time evolution of the momentum distribution of non-equilibrated heavy quarks due to the interaction with the equilibrated QGP during the time interval τ i < τ < τ HQ can be calculated by the Fokker-Planck equation. Evolution of the expanding QGP is described by relativistic hydrodynamics [32,33].
To derive the Fokker-planck equation, let us start from the Boltzmann transport equation, which in covariant form is [34]: where p µ is the four-momentum of the heavy quarks (HQs), C{f } is the collision term and f (x, p) is the momentum distribution function of the heavy quarks. For uniform plasma, f will lose its dependence on x, so the Boltzmann equation becomes To simplify the non-linear integro-differential Boltzmann equation, we employ the Landau approximation. This means that we consider soft scatterings with small momentum transfer compared to the particle momentum p, and neglect the interaction of heavy quarks with each other. According to these assumptions the evolution of momentum distribution of heavy quarks is given by where'd, A i (p) and B ij (p) are drag and diffusion coefficients respectively, defined as: Here, we considered the elastic scattering of the HQs with gluons, light quarks, and corresponding anti-quarks. All these interactions contribute to determine the collision rate w(p, k) [35]. For an isotropic QGP one can write A i = p i A(p) and B i,j = D(p)δ i,j . In the simplest form, we assume that the drag and diffusion coefficients are slowly varying functions of momentum, so Equation (3) reduces to the Fokker-Planck equation: The drag coefficient is an important quantity carrying information about the dynamics of elastic collisions. It is expected that the drag coefficient should be determined by the properties of the bath and not by the nature of the HQ [36]. After heavy ion collision, the charm quarks are produced at a time scale of about 0.07f m/c. This time is about 0.02f m/c for bottom quarks. These heavy quarks will propagate through the deconfined matter (QGP) while losing energy.
In this work we use the drag coefficient as follows: If coupling between HQs and background QGP is weak, then the diffusion coefficient can be calculated using the Einstein relation, D(p) = M T A(p), where M is the mass of the heavy quark and T is the thermal bath temperature [37,38]. The first perturbative estimation of the energy loss in a QGP was proposed by Bjorken for heavy quarks [39], which is an important quantity. It may be noted that the jet quenching, drag force, damping rate of particles in the plasma and etc., are determined directly in terms of the energy loss. Furthermore, since the calculated energy loss is proportional to the heavy fermion mass in quark-gluon plasma, it is used as an important tool to identify particles and give more information about the properties of QGP. The collisional energy loss of a heavy quark with energy E and mass M in a thermalized medium with temperature T has been calculated in several papers with different approaches [40]. In our calculations, the collisional energy loss of HQs in QGP by considering "Hard and Soft Thermal Loops" are given as follows [41].
where v is the HQ speed, α s is the strong coupling constant, n f is the number of quark flavors in the medium, E and M are energy and mass of the propagated HQ respectively, m g = (1 + n f /6)g 2 T 2 /3 is the thermal gluon mass, g = √ 4πα s is the gauge coupling parameter and B(v) is a smooth velocity function, which can be taken approximately as 0.7 [41].
The collisional energy loss of a HQ in a thermalized QGP is calculated by "Hard Thermal Loop" approximation as follows [42]: where m D = 8αs π (3 + n f )T 2 is the Debye mass. Finally the total energy loss due to both collision and radiation precesses is [43]: where ξ = 9E 2π 3 T . Using equation (8) − (11), we will calculate the drag and diffusion coefficients for HQ energy loss during the interaction with QGP medium. In the next section we will find the solution of the Fokker-plank equation, eq(6).

III. The Iterative Laplace Transform Method
In this section we want to introduce a new method for solving the FP equation, namely the "Iterative Laplace Transform Method" (ILTM) [44]. Many problems in physics and engineering can be successfully modelled by fractional differential equations (FDE). The FP equation indeed, is an important example of FDEs. Standard methods like finitedifference method are not applicable for such differential equations with non-zero values on the boundaries of initial conditions. Hence should we be looking for an effective way to solve FDEs. This new method was proposed by Daftardar-Gejji and Jafari to seek numerical solutions of nonlinear functional equations. This method is called the Iterative Laplace Transform Method, which is a combination of two powerful methods, namely, the Laplace transform method and the Iterative method. The method gives numerical solutions in the form of convergent series with easily computable components. The most outstanding feature of this method is that it provides an analytical solution using the initial conditions only, without any discretization or restrictive assumptions.
In general, the Fokker-Planck equation can be written as [44]: Here D α t (.), D β p (.), D 2β p (.) are the Caputo fractional derivative with respect to t and p. In our calculations we use the following definitions: I. The Caputo fractional derivative of function (p, t) is defined as II. The Laplace transform of f (t) is defined as: III. Laplace transform of D α t f (p, t) is given by where f (k) (p, 0) is the kth derivative of f (p, t) at t = 0.
To illustrate how the iterative Laplace transform method works, the general space-time fractional partial differential equation is considered where and f (p, t) will be determined later.
First of all we take the Laplace transform of both sides of (16) After simplification, we have By inverse Laplace transforms of both sides of Eq. (19) we have: and The method gives numerical solutions in the form of convergent series and By substituting (22) and (23) into (20) we have Then we set Therefore the solution of (16)with initial condition (17)is given by: If we set α = β = 1 in Equation (26) we will reach the FPE. The drag and diffusion coefficients are explicitly timely independent parameters. However we recalculate these parameters before every time steps of our simulation. In this situation we have: where We solved the FPE using both finite difference method and ILTM with different types of initial conditions and compared their results. For initial conditions with non vanishing values at the borderer, the ILTM successfully solves the FPE while results of finite difference method diverges and therefore it failes to find the solution. It may be noted that momentum distribution functions of HQs, generally find their maximum at (p = 0). We have taken momentum distribution functions of charm and bottom quarks at √ s = 200GeV in Au-Au and also P-P collisions from "The Durham HepData Project" database [45] which are simulated data. ILTM has been performed up to m = 2, which means we calculate f 0 , f 1 and f 2 in (28).

IV. Time evolution of HQ distribution functions in QGP
The dissipation coefficients of the thermal bath and initial momentum distribution of HQs, are basic inputs required to solve the FPE. Time evolution of dynamical parameters (like temperature, viscosity and so on) should be taken from suitable models. In our calculation the time dependence of temperature is considered as [46] T (t) = τ where τ 0 and T 0 are the initial time and initial temperature that the background of the partonic system has attained in local equilibrium. We have taken τ 0 = 0.33 fm/c and T 0 = 0.373 GeV in our simulations [46]. The viscosity effects have not considered for simplification of our problem. Here we considered a fixed value for coupling constant as α s = 0.3. It is clear that the background is not stationary, it expands and cools down with time, so for more accurate results one can consider a strong coupling, α s , running with HQ mass and/or the medium temperature. For calculating the time evolution of the HQ momentum distribution by solving the FPE using ILTM, we need three inputs: I) initial distribution function, II) drag coefficient A(p) which can be calculated by Inserting (11) in the Equation (7) and III) the diffusion coefficient as D = A(p)T M . Finally, we find the momentum distribution of HQs at transition temperature T c = 175 MeV. It is clear that more advanced relations and more suitable assumptions can produce better results which can be found in further works. Figure 1 demonstrates the drag coefficient for b and c quarks due to collision and radiation separately. This figure clearly shows that the energy loss for c quark is greater than that for b quark. Figure 2 presents the HQ momentum distributions at initial time (τ = 0.33 fm/c) and at final time, when the QGP temperature becomes T = 0.175 GeV. A general comparison of momentum distributions before and after the evolution indicates that the probability of finding quarks with lower momentums after the evolution is greater than the probability of finding low momentum quarks at initial time. It may be noted that momentum distributions can not be compared directly, because they describe quark distributions at different values of total energy. This figure also shows that energy loss due to suppression of c quark is greater than that for b quark.
V. Calculating the nuclear suppression factor R AA (P T ) One of the most important experimental observables to quantify the nuclear medium is the depletion of high p T particles (D and B mesons or non-photonic single e − spectra resulting from semi leptonic decays of hadrons containing charm and bottom quarks) produced in Nucleus-Nucleus collisions with respect to those produced in proton-proton collisions which is formulated in nuclear suppression factor R AA . The nuclear suppression factor is calculated by solving relativistic hydrodynamic equations for the background medium in QGP phase, simultaneously with the FP equation for the HQs.The properties of the background bath affect the energy loss of HQs through the drag and diffusion coefficients. The energy loss of HQs is sensitive to initial conditions and the model used to describe the dynamics of the system. Here the nuclear modification factor (R AA ) in relativistic heavy ion collisions is calculated for demonstrating the power of ILTM applied to solve the FP equation.
We now proceed with the presentation of our numerical simulation with the formalism described and the inputs mentioned in the previous section. Then we present the comparison between the results obtained with solving FPE by ILTM approach by the experimental data. There are several free parameters in the problem, which should be initiated or chosen to find an acceptable fitting the results on the reference data. Equilibration time after the collision, equilibration temperature, contribution of collisional and radiation dissipations on drag and diffusion coefficients, effective light quark flavors n f , the charm and the bottom mass, hadron multiplicity at mid-rapidity dN dη and final temperature are some of needed parameters.
To make a realistic connection between simulation results and the experimental data, we have to implement the hadronization mechanisms at the quark-hadron transition temperature in order to find the non-photonic single electron spectra originating from the decays of heavy flavoured mesons (D and B). The solution of the FP equation (described in the previous section) have been convoluted with the fragmentation functions of the heavy quarks to obtain the P T distribution of the heavy mesons. We have used the following heavy quark fragmentation [46,47]: We have assumed that the electron spectrum produced through semileptonic decay of D and B-mesons is proportional to the heavy meson spectra. Also we have not considered the expansion of the QGP bath. These assumptions create some expected deviation from the experimental data.
Similar simulations also have been performed to find the electron spectrum created in the P-P collisions. For this purpose HQ distribution functions of proton at √ s = 200GeV have been taken as initial conditions to solve FPE again. The ratio of these two quantities is proportional to the nuclear suppression factor, R AA measured in experiments: A small value of R AA indicates a strong suppression and therefore, large energy loss of heavy quarks in the medium. It is clear that this ratio is equal to one in the absence of any medium. Figure 3 presents contribution of collision and radiation energy loss in suppression factor R AA separately. Displayed experimental data has been obtained from the PHENIX collaborations for Au + Au collisions at √ s N N = 200GeV . This figure shows that suppression at higher values of P T can not be explained by collisional energy loss. In other words, contribution of radiation in energy loss at higher P T is more important than collisional energy loss only. Figure 4 demonstrates our best simulation result by considering both collision and radiation energy loss in calculating the drag and diffusion coefficient to solve the FPE. The drag coefficient has been considered as A total = A collison + 1.7A radiation . This figure shows a very good agreement between results and experimental data at least for P T < 5 GeV. It may be noted that the final temperature has been taken T c = 0.16 GeV and an offset 0.05 has been added too. As mentioned before we have not considered some details of the real situation in our simulations.

IV. Conclusions and remarks
The "Iterative Laplace Transform Method" has been introduced as an effective method to solve the Fokker-Planck equation for initial conditions with non-zero values at the boundaries. A numerical algorithm for solving FPE is presented in this paper. Our calculations show that this method is able to solve time evolution of distribution function of HQs successfully. To verify the ability of ILTM and demonstrating its application, we calculated the nuclear suppression factor R AA . It is shown that there is a good agreement between simulation results and reported experimental data. Therefore many problems can be investigated using the ILTM. Effects of a running coupling constant, HQ mass, dead cone effect and LPM effects, QGP viscosity, different equations of state as models of QGP are subjects which can be studied in further works.