Heavy Quarks in Turbulent QCD Plasmas

TL;DR

This paper develops a Fokker-Planck model to study how heavy quarks interact with turbulent quark-gluon plasmas, revealing differences in momentum broadening and energy loss compared to equilibrium conditions.

Contribution

It introduces a novel approach to model heavy quark dynamics in turbulent QCD plasmas using classical fields, extending beyond equilibrium assumptions.

Findings

Heavy quarks experience different momentum broadening in turbulent vs. equilibrium plasmas.

Energy loss of heavy quarks is affected by the presence of long wavelength chromodynamic fields.

Turbulent plasmas modify heavy quark transport properties compared to equilibrated plasmas.

Abstract

The quark-gluon plasma, which is produced at an early stage of ultrarelativistic heavy-ion collisions, is expected to be initially strongly populated with chromodynamic fields. We address the question how heavy quarks interact with such a turbulent plasma in comparison with an equilibrated one of the same energy density. For this purpose we derive a Fokker-Planck transport equation of heavy quarks embedded in a plasma of light quarks and gluons. We first discuss the equilibrium plasma and then the turbulent one applying the same approach, where the heavy quarks interact not with the plasma constituents but rather with the long wavelength classical fields. We first consider the three schematic models of isotropic trubulent plasma and then the simplified model of glasma with the chromodynamic fields only along the beam direction. The momentum broadening and collisional energy loss of a test heavy quark are computed and compared to those of equilibrium plasma of the same energy density.

Tables

Table 1. Table 1: The coefficients X L ​ ( v ) subscript 𝑋 𝐿 𝑣 X_{L}(v) and X T ​ ( v ) subscript 𝑋 𝑇 𝑣 X_{T}(v) in the models of turbulent plasmas

		$\frac{1}{g^2{C}_F}{X}_L\left(v\right)$		$\frac{1}{g^2{C}_F}{X}_T\left(v\right)$
Gauss E & B		$\sqrt{\frac{\pi }{2}}\frac{M_E{\sigma }_t{\sigma }_r}{\sqrt{{\sigma }_r^2+v^2{\sigma }_t^2}}$		$\sqrt{\frac{\pi }{2}}\frac{\left(M_E+v^2{M}_B\right){\sigma }_t{\sigma }_r}{\sqrt{{\sigma }_r^2+v^2{\sigma }_t^2}}$
Gauss A		$\sqrt{\frac{\pi }{2}}\frac{M_A{\sigma }_t{\sigma }_r{v}^2}{{\left({\sigma }_r^2+v^2{\sigma }_t^2\right)}^{3/2}}$		$\sqrt{\frac{\pi }{8}}\frac{M_A{\sigma }_t{v}^2}{{\sigma }_r}\left[\frac{3}{{\left({\sigma }_r^2+v^2{\sigma }_t^2\right)}^{1/2}}-\frac{{\sigma }_r^2+3{\sigma }_t^2{v}^2}{{\left({\sigma }_r^2+v^2{\sigma }_t^2\right)}^{3/2}}\right]$
Stationary A		0		$\frac{ℳv}{8\pi }\left[\frac{1}{2}{k}_{\max }^2-\frac{1}{2}{\mu }^2\ln \left(\frac{k_{\max }^2+{\mu }^2}{{\mu }^2}\right)\right]$

Table 2. Table 2: The transport coefficients d ​ E d ​ x 𝑑 𝐸 𝑑 𝑥 \frac{dE}{dx} and q ^ ^ 𝑞 \hat{q} in the models of turbulent plasmas.

		$-\frac{1}{g^2{C}_F}\frac{dE}{dx}$		$\frac{1}{g^2{C}_F}\hat{q}$
Gauss E & B		$\sqrt{\frac{\pi }{2}}\frac{M_E{\sigma }_t{\sigma }_{r}v}{{\left({\sigma }_r^2+v^2{\sigma }_t^2\right)}^{1/2}T}$		$\sqrt{8\pi }\frac{\left(M_E+v^2{M}_B\right){\sigma }_t{\sigma }_r}{v\sqrt{{\sigma }_r^2+v^2{\sigma }_t^2}}$
Gauss A		$\sqrt{\frac{\pi }{2}}\frac{M_A{\sigma }_t{\sigma }_r{v}^3}{{\left({\sigma }_r^2+v^2{\sigma }_t^2\right)}^{3/2}T}$		$\sqrt{\frac{\pi }{2}}\frac{M_A{\sigma }_{t}v}{{\sigma }_r}\left[\frac{3}{{\left({\sigma }_r^2+v^2{\sigma }_t^2\right)}^{1/2}}-\frac{{\sigma }_r^2+3v^2{\sigma }_t^2}{{\left({\sigma }_r^2+v^2{\sigma }_t^2\right)}^{3/2}}\right]$
Stationary A		0		$\frac{ℳ}{2\pi }\left[\frac{1}{2}{k}_{\max }^2-\frac{1}{2}{m}^2\ln \left(\frac{k_{\max }^2+m^2}{m^2}\right)\right]$

Table 3. Table 3: Parametric estimates of the coefficients X T ​ ( v ) subscript 𝑋 𝑇 𝑣 X_{T}(v) in the four plasma models under consideration. The third column takes into account that σ t − 1 ∼ k max ∼ g ​ T similar-to superscript subscript 𝜎 𝑡 1 subscript 𝑘 max similar-to 𝑔 𝑇 \sigma_{t}^{-1}\sim k_{\rm max}\sim gT .

		$X_T\left(v\right)$		$X_T\left(v\right)$
equilibrium		$g^4{T}^3$		$g^4{T}^3$
Gauss E & B		$g^2\tau T^4$		$g{T}^3$
Gauss A		$g^2\tau T^4{v}^2$		$g{T}^3{v}^2$
Stationary A		$g^2\frac{T^4}{k_{\max }}v$		$g{T}^{3}v$

Equations

$$\big{(}D^{0}+{\bf v}\cdot{\bf D}\big{)}Q(t,{\bf r},{\bf p})-{1\over 2}\big{\{}{\bf F}(t,{\bf r}),\nabla_{p}Q(t,{\bf r},{\bf p})\big{\}}=0,$$

$$Q(t,\mathbf{r},\mathbf{p})=⟨Q(t,\mathbf{r},\mathbf{p})⟩+\delta Q(t,\mathbf{r},\mathbf{p}),$$

$$⟨Q(t,\mathbf{r},\mathbf{p})⟩=n(t,\mathbf{r},\mathbf{p})\text{mathds}1,$$

$$|⟨Q⟩|\gg |\delta Q|,|{\nabla }_p⟨Q⟩|\gg |{\nabla }_p\delta Q|,$$

$$\left|\frac{\partial \delta Q}{\partial t}\right|\gg \left|\frac{\partial ⟨Q⟩}{\partial t}\right|,|\nabla \delta Q|\gg |\nabla ⟨Q⟩|.$$

$$⟨\mathbf{E}(t,\mathbf{r})⟩=⟨\mathbf{B}(t,\mathbf{r})⟩=0.$$

$$\mathcal{D}\delta Q(t,\mathbf{r},\mathbf{p})-\mathbf{F}(t,\mathbf{r})\cdot {\nabla }_{p}n(t,\mathbf{r},\mathbf{p})=0,$$

$${\cal D}\,n(t,{\bf r},{\bf p})-\frac{1}{N_{c}}{\rm Tr}\big{\langle}{\bf F}(t,{\bf r})\cdot\nabla_{p}\delta Q(t,{\bf r},{\bf p})\big{\rangle}=0.$$

$$\delta Q(t=0,\mathbf{r},\mathbf{p})=\delta Q_0(\mathbf{r},\mathbf{p}),$$

$$f(\omega ,\mathbf{k})={\int }_0^{\infty }dt\int d^{3}r{e}^{i(\omega t-\mathbf{k}\cdot \mathbf{r})}f(t,\mathbf{r}).$$

$$f(t,\mathbf{r})={\int }_{-\infty +i\sigma }^{\infty +i\sigma }\frac{d\omega }{2\pi }\int \frac{d^{3}k}{(2\pi {)}^3}{e}^{-i(\omega t-\mathbf{k}\cdot \mathbf{r})}f(\omega ,\mathbf{k}),$$

$$\delta Q(\omega ,\mathbf{k},\mathbf{p})=i\frac{\mathbf{F}(\omega ,\mathbf{k})\cdot {\nabla }_{p}n(\mathbf{p})+\delta Q_0(\mathbf{k},\mathbf{p})}{\omega -\mathbf{k}\cdot \mathbf{v}},$$

$$\delta Q(t,{\bf r},{\bf p})=\int_{0}^{t}dt^{\prime}\>{\bf F}\big{(}t^{\prime},{\bf r}-{\bf v}(t-t^{\prime})\big{)}\cdot\nabla_{p}n({\bf p})+\delta Q_{0}({\bf r}-{\bf v}t,{\bf p}),$$

$$\displaystyle\big{\langle}{\bf F}(t,{\bf r})\cdot\nabla_{p}\delta Q(t,{\bf r},{\bf p})\big{\rangle}=\int_{0}^{t}dt^{\prime}\>\nabla_{p}^{i}\big{\langle}F^{i}(t,{\bf r})\,F^{j}\big{(}t^{\prime},{\bf r}-{\bf v}(t-t^{\prime})\big{)}\big{\rangle}\nabla_{p}^{j}n({\bf p})+\nabla_{p}^{i}\big{\langle}F^{i}(t,{\bf r})\delta Q_{0}({\bf r}-{\bf v}t,{\bf p})\big{\rangle}.$$

$$\frac{1}{N_{c}}{\rm Tr}\big{[}\big{\langle}F^{i}(t,{\bf r})\delta Q_{0}({\bf r}-{\bf v}t,{\bf p})\big{\rangle}\big{]}=Y^{i}({\bf v})\,n({\bf p}),$$

$$X^{ij}({\bf v})\equiv\frac{1}{N_{c}}\int_{0}^{t}dt^{\prime}\>{\rm Tr}\big{[}\big{\langle}F^{i}(t,{\bf r})F^{j}\big{(}t^{\prime},{\bf r}-{\bf v}(t-t^{\prime})\big{)}\big{\rangle}\big{]},$$

$$\Big{(}{\cal D}-\nabla_{p}^{i}X^{ij}({\bf v})\nabla_{p}^{j}-\nabla_{p}^{i}Y^{i}({\bf v})\Big{)}n(t,{\bf r},{\bf p})=0.$$

$$⟨E_a^i(t_1,{\mathbf{r}}_1)E_a^j(t_2,{\mathbf{r}}_2)⟩⟶⟨E_a^i(t_1,{\mathbf{r}}_1){\Omega }_{ab}(t_1,{\mathbf{r}}_1|t_2,{\mathbf{r}}_2)E_b^j(t_2,{\mathbf{r}}_2)⟩,$$

$$\Omega(t_{1},{\bf r}_{1}|t_{2},{\bf r}_{2})={\cal P}\exp\Big{[}ig\int_{(t_{2},{\bf r}_{2})}^{(t_{1},{\bf r}_{1})}ds_{\mu}A^{\mu}_{c}(s)T^{c}\Big{]}.$$

$$\Omega (t_1,{\mathbf{r}}_1|t_2,{\mathbf{r}}_2)⟶\mathcal{U}(t_1,{\mathbf{r}}_1)\Omega (t_1,{\mathbf{r}}_1|t_2,{\mathbf{r}}_2){\mathcal{U}}^T(t_2,{\mathbf{r}}_2),$$

$$X^{ij}(\mathbf{v})$$

$$Y^i(\mathbf{v})$$

$$\displaystyle X_{L}(v)=\frac{v^{i}v^{j}}{{\bf v}^{2}}X^{ij}({\bf v}),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}X_{T}(v)=\frac{1}{2}\Big{(}\delta^{ij}-\frac{v^{i}v^{j}}{{\bf v}^{2}}\Big{)}X^{ij}({\bf v}),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}Y(v)=\frac{v^{i}}{{\bf v}^{2}}Y^{i}({\bf v}).$$

$$n^{\rm eq}({\bf p})\sim\exp\Big{(}-\frac{E_{\bf p}}{T}\Big{)},$$

$$X^{ij}(\mathbf{v})\frac{v^j}{T}=Y^i(\mathbf{v}),$$

$$X_L(v)\frac{1}{T}=Y(v).$$

$$\bigg{(}{\cal D}-X\Big{(}\nabla_{p}^{2}+\frac{1}{T}\nabla_{p}\cdot{\bf v}\Big{)}\bigg{)}n(t,{\bf r},{\bf p})=0,$$

$$\frac{⟨\Delta p^i⟩}{\Delta t}$$

$$\frac{⟨\Delta p^i\Delta p^j⟩}{\Delta t}$$

$$\displaystyle\frac{\langle\Delta E\rangle}{\Delta t}=\frac{1}{\Delta t}\Big{\langle}\frac{{\bf p}\cdot\Delta{\bf p}}{\sqrt{m^{2}+{\bf p}^{2}}}\Big{\rangle}=v^{i}\frac{\langle\Delta p^{i}\rangle}{\Delta t}.$$

$$\frac{dE}{dx}=\frac{v^i}{v}\frac{⟨\Delta p^i⟩}{\Delta t}=-\frac{v^i}{v}{Y}^i(\mathbf{v}),$$

Full text

I Introduction

The early stage of relativistic heavy-ion collisions is least known because there are hardly any experimentally accessible signals of the phase. Nevertheless one expects that the quark-gluon plasma, which is produced in the collisions, is initially strongly populated with chromodynamic fields. Within the framework of the Color Glass Condensate (CGC) approach, see e.g. the review Gelis:2012ri , color charges of partons confined in the colliding nuclei act as sources of long wavelength chromodynamic fields which can be treated classically because of large occupation numbers of the soft modes. Since the density of color charges per transverse area of heavy nuclei is large, the corresponding momentum scale $Q_{s}$ is expected to be significantly bigger than the QCD scale parameter $\Lambda_{\rm QCD}$. Consequently, the coupling constant $\alpha_{s}$ is presumably sufficiently small and perturbative methods are applicable. The system, however, is rather strongly interacting because of the high-amplitude fields present in the system.

A momentum anisotropy of the early stage quark-gluon plasma makes it unstable with respect to chromomagnetic modes which in turn cause a spontaneous generation of the fields, as explained at length in the review article Mrowczynski:2016etf . Therefore, the effect of strong fields is further enhanced. Following the terminology of electromagnetic plasma, we call such a nonequilibrium system of fields as the turbulent plasma meaning that numerous modes are excited in the system. In the CGC approach the non-equilibrium system of fields from the early stage of relativistic heavy-ion collisions is called glasma Gelis:2012ri and it can be treated as a specific realization of the turbulent QCD plasma. Leaving aside the mechanism of field generation and its detailed structure, one asks what are the transport properties of turbulent plasmas. We are specifically interested how heavy quarks - charm or beauty - behave in such a system when compared to the equilibrium plasma of the same energy density.

Heavy quarks are often treated as a probe of strongly interacting matter created in relativistic heavy-ion collisions, see e.g. the review Prino:2016cni . Thanks to their large masses the quarks are produced only at the earliest stage of the collision due to hard interactions of partons from incoming nuclei. Later on they propagate through a surrounding medium testing the entire history of the system. It has been long believed that the interaction of heavy quarks is significantly weaker than that of light quarks or gluons but experimental data clearly contradict the expectation. As discussed in the review Prino:2016cni , the behavior of mesons containing a heavy quark is rather similar to that of light mesons at both small and large transverse momenta. The problem is not fully resolved.

The medium created in relativistic heavy-ion collisions evolves fast towards the locally equilibrated quark-gluon plasma which expands hydrodynamically and ultimately is converted into a hadron gas. Final momentum spectra of heavy quarks are mostly shaped in the long-lasting equilibrium phase which is relatively well understood. An effect of a pre-equilibrium phase is often entirely ignored but this transient phase can significantly influence heavy-quark spectra because of its high density. Non-equilibrium calculations recently performed in a framework of kinetic theory Das:2017dsh confirm the suggestion. However, we are interested even in the earlier phase when the medium is not described in terms of quasi-particles, as in a kinetic theory, but rather as a system dominated by classical fields which is the turbulent plasma.

A simple parametric estimate suggests that the interaction of heavy quarks in a turbulent plasma is much stronger than in the equilibrium one, if the coupling constant $g$ is small. The momentum broadening parameter $\hat{q}$, for example, is of order $g^{4}$ in equilibrium plasmas. Since the quark of interest actually interacts with soft gluons emitted by plasma constituents, one can think that the factor $g^{4}$ is composed of two pieces of $g^{2}$. The first one is related to the gluon emission and the second one to the gluon absorption. If soft chromodynamic fields are present in the plasma, the interaction of the quark should be rather of order $g^{2}$ than $g^{4}$. In Sec. VI we argue that $\hat{q}$ is indeed not of the order $g^{4}$, not even $g^{2}$ but presumably of the order $g$ in a turbulent plasma.

Because of their big masses, relaxation times of heavy quarks, which are produced in relativistic heavy-ion collisions, are expected to be significantly longer than that of light quarks and gluons. When an equilibrium or, more generally, a stationary state is reached by light quarks and gluons, heavy quarks need some extra time to adjust to the state of the plasma. Such a situation is naturally described in terms of the Fokker-Planck transport equation which was indeed repeatedly applied to heavy quarks in Moore:2004tg ; Svetitsky:1987gq ; vanHees:2004gq ; Mustafa:2004dr . The equation is usually derived from the Boltzmann equation by applying the so-called diffusion approximation to the collision term LP81 . The approximation assumes that the momentum transfer to the heavy quark in every collision is much smaller than the quark momentum.

The aim of this paper is threefold. In the first part we rederive the Fokker-Planck equation of heavy quarks which do not interact with plasma constituents but rather with soft classical fields present in the plasma. Specifically, we apply the so-called quasi-linear theory known from the electromagnetic plasma Ved61 ; LP81 . The theory assumes that the distribution function can be decomposed into a large but slowly varying regular part and a small fluctuating or turbulent one which oscillates fast. The average over a statistical ensemble of the turbulent part is assumed to vanish and thus the average of the distribution function equals its regular part. The turbulent contribution to the distribution function obeys the collisionless transport equation while the transport equation of the regular part is determined by the fluctuation spectra which provide the collision term. The derivation presented here closely follows the procedure which was developed for QCD in Mrowczynski:2009gf , where, however, only the longitudinal chromoelectric field was taken into account and here the complete chromodynamic field is considered. The equilibrium correlation functions of chromodynamic fields, which are needed to obtain the quasi-linear transport equations, were derived in Mrowczynski:2008ae .

Our second aim is to confront the equilibrium plasma with the turbulent one. Therefore, we consider three models of isotropic turbulent plasma in the second part of the paper. Postulating a form of the correlation functions of chromodynamic fields, we derive the coefficients of the Fokker-Planck equation which can be related to the energy loss, momentum broadening and diffusion coefficient of heavy quarks in the plasma. The transport coefficients of turbulent plasma are compared to those of the equilibrium one at the same energy density.

The third aim is to study an evolution of heavy quarks at the earliest stage of relativistic have-ion collisions. Since the chromoelectric and chromomagnetic fields spanned between the receding nuclei are initially mostly parallel to the beam direction, we model the glasma with the boost invariant correlation functions of longitudinal fields. The energy loss and momentum broadening of heavy quarks are computed, assuming that all energy of the glasma is accumulated in the longitudinal fields.

At the end of the introductory remarks we note that an approach similar to ours, which was also inspired by the electromagnetic plasma studies Dupree:1966 , was formulated in Asakawa:2006jn , see Majumder:2007zh ; Chandra:2015gma as well. We also mention an attempt Song:2015jmn to study transport of heavy quarks in a plasma populated by strong chromodynamic fields. Unfortunately, the paper is flawed as the framework of an isotropic Langevin approach is applied to anisotropic plasmas.

Throughout the paper we use the natural system of units with $c=\hbar=k_{B}=1$; our choice of the signature of the metric tensor is $(+---)$. Lorentz indices are denoted with $\mu,\,\nu=0,\,1,\,2,\,3$ and $i,\,j=1,\,2,\,3$ label the Cartesian coordinates $x,\,y,\,z$. The color indices of the adjoint representation of ${\rm SU}(N_{c})$ gauge group are $a,\,b=1,\,2,\dots N_{c}^{2}-1$.

II Derivation of Fokker-Planck equation

Our derivation of the Fokker-Planck equation of heavy quarks embedded in quark-gluon plasma starts with the transport equation of the Vlasov form

[TABLE]

where the distribution function $Q(t,{\bf r},{\bf p})$ of heavy quarks is the $N_{c}\times N_{c}$ hermitian matrix which belongs to the fundamental representation of the SU($N_{c}$) group. The distribution function depends on the time ($t$), position (${\bf r}$) and momentum (${\bf p}$) variables. There is no explicit dependence on the time-like component of the four-momentum $p^{\mu}=(p^{0},{\bf p})$ as the distribution function is assumed to be non-zero only for momenta obeying the mass-shell constraint that is $p_{0}=E_{\bf p}=\sqrt{{\bf p}^{2}+m^{2}}$. The quark velocity equals ${\bf v}={\bf p}/E_{\bf p}$ and $D^{\mu}\equiv(D^{0},{\bf D})\equiv\partial^{\mu}-ig[A^{\mu}(x),\cdots\;]$ with $A^{\mu}(x)$ being the chromodynamic potential in the fundamental representation. The mean-field term of the transport equation (1) is expressed through the color Lorentz force {\bf F}(t,{\bf r})\equiv g\big{(}{\bf E}(t,{\bf r})+{\bf v}\times{\bf B}(t,{\bf r})\big{)} with the chromoelectric ${\bf E}(t,{\bf r})$ and chromomagnetic ${\bf B}(t,{\bf r})$ fields also belonging to the fundamental representation. The symbol $\{\dots,\dots\}$ denotes the anticommutator. The derivation of the transport equation (1) is discussed in detail in the review Mrowczynski:2016etf .

Further on we assume that the chromodynamic fields and the distribution function which enter the transport equation (1) can be decomposed into a regular and fluctuating or turbulent component. The distribution function is thus written down as

[TABLE]

where $\langle\cdots\rangle$ denotes ensemble average; $\langle Q(t,{\bf r},{\bf p})\rangle$ is called the regular part while $\delta Q(t,{\bf r},{\bf p})$ is called the fluctuating or turbulent one. It directly follows from Eq. (2) that $\langle\delta Q\rangle=0$. The regular contribution is assumed to be color neutral or white, and it is expressed as

[TABLE]

where $\mathds{1}$ is the unit matrix in color space. Since the distribution function transforms under gauge transformations as $Q\rightarrow U\,QU^{\dagger}$, where $U$ is the transformation matrix, the regular contribution of the form (3) is gauge independent. We also assume that

[TABLE]

but at the same time

[TABLE]

What concerns the chromodynamic fields, we assume in accordance with Eq. (3) that their regular parts vanish and thus

[TABLE]

We substitute the distribution function (2) into the transport equation (1) and linearize the equations in the fluctuating contributions. Thus we get the equation

[TABLE]

where ${\cal D}\equiv\frac{\partial}{\partial t}+{\bf v}\cdot\nabla$ is the substantial or material derivative.

Now we substitute the distribution functions (2) into the transport equations (1) but instead of linearizing the equation in the fluctuating contributions, we take the ensemble average of the resulting equation and trace over the color indices. Thus we get

[TABLE]

Since the regular part of distribution function is assumed to be color neutral, see Eq. (3), the term ${\rm Tr}[\langle{\bf F}\cdot\nabla_{p}n\rangle]$ vanishes because the fields ${\bf E},\;{\bf B}$ are traceless. The trace over color indices also cancels the terms originating from covariant derivatives like ${\rm Tr}\langle[A^{\mu},\delta Q]\rangle$. We finally note that the trace ${\rm Tr}[\langle{\bf F}\cdot\nabla_{p}\delta Q\rangle]$ is gauge independent as the regular distribution function $n(t,{\bf r},{\bf p})$ is.

Now, we are going to write down the transport equation (8) in the Fokker-Planck form. For this purpose we observe that due to the condition (5), the space-time dependence of the regular distribution function can be neglected in the linearized transport equation (7) and then, the equation becomes easily solvable. We solve it with the initial condition

[TABLE]

using the one-sided Fourier transformation defined as

[TABLE]

The inverse transformation is

[TABLE]

where the real parameter $\sigma>0$ is chosen in such a way that the integral over $\omega$ is taken along a straight line in the complex $\omega-$plane, parallel to the real axis, above all singularities of $f(\omega,{\bf k})$.

The linearized transport equation (7), which is converted into the algebraic equation by means of the one-sided Fourier transformation, is solved as

[TABLE]

We stress that although we have ignored the (weak) frequency and wave number dependence of the regular distribution $n$, the fields ${\bf E}(\omega,{\bf k}),\;{\bf B}(\omega,{\bf k})$ retain their full frequency and wave number dependence in the expression (12). Inverting the one-sided Fourier transformation, one finds the solution of the linearized transport equation as

[TABLE]

where we assumed that ${\bf E}(\omega,{\bf k})$ and ${\bf B}(\omega,{\bf k})$ are analytic functions of $\omega$.

With the help of the solution (13), the force term in the transport equation (8) becomes

[TABLE]

The second term in the r.h.s. of Eq. (14) can be manipulated to the form

[TABLE]

which is effectively the definition of the vector $Y^{i}({\bf v})$. We also introduce the tensor

[TABLE]

and we note that, as explained in the subsequent sections, $X$ and $Y$ become time independent for a sufficiently long $t$. Then, the transport equation (8) can be written as the Fokker-Planck equation

[TABLE]

Since the distribution function $n(t,{\bf r},{\bf p})$ carries no information about color degrees of freedom, the function is gauge invariant, and consequently $X^{ij}({\bf v})$ and $Y^{i}({\bf v})$ should be gauge invariant as well. However, one observes that $X^{ij}({\bf v})$ and $Y^{i}({\bf v})$ as defined by Eqs. (15) and (16) are gauge dependent because the traces are of nonlocal quantities in the definitions (15) and (16). The starting transport equation (1) is gauge covariant but the linearization procedure breaks the covariance because the covariant derivative is replaced by the normal one. Consequently, the solution (13) is not gauge covariant – the right-hand side of Eq. (13) transforms differently under local gauge transformations than the left-hand side. To cure the problem, one modifies the solution (13) by means of the link operator which is also called the gauge parallel transporter, see e.g. Sec. IIIE of the review article Mrowczynski:2016etf . Then, the modified solution obeys Eq. (8) with the covariant derivative instead of the normal one. Let us briefly discuss the procedure in a context of the Fokker-Planck equation (17).

According to Eq. (16), the tensor $X^{ij}({\bf v})$ is determined by the traces of the field correlation functions like $\langle E^{i}_{a}(t_{1},{\bf r}_{1})\,E^{j}_{a}(t_{2},{\bf r}_{2})\rangle$ where chromodynamic fields are written in the adjoint representation of the ${\rm SU}(N_{c})$ group which is used further on. The trace becomes gauge invariant under the replacement

[TABLE]

where $\Omega_{ab}(t_{1},{\bf r}_{1}|t_{2},{\bf r}_{2})$ is the link operator defined as

[TABLE]

Here $T^{c}$ is the adjoint representation generator of the ${\rm SU}(N_{c})$ group and ${\cal P}$ denotes the ordering along the path connecting the points $(t_{2},{\bf r}_{2})$ and $(t_{1},{\bf r}_{1})$. Since the fields transform as vectors under the local gauge transformation ${\cal U}(t,{\bf r})$ and the link transforms as

[TABLE]

one checks that the trace of the correlation function which includes the link is indeed gauge invariant. Consequently, the tensor $X^{ij}({\bf v})$ is gauge invariant. Analogously one achieves the gauge invariance of the vector $Y^{i}({\bf v})$. Further on, whenever any nonlocal correlation function shows up a presence of the link operator is implicitly assumed even so it is not explicitly written.

III Fokker-Planck equation

Although this is a textbook material we briefly discuss here the Fokker-Planck equation (17). We first note that in the isotropic plasma the tensor $X^{ij}({\bf v})$ and vector $Y^{i}({\bf v})$ both depend on a single vector that is the heavy-quark velocity ${\bf v}$. Therefore, they can be written as

[TABLE]

where $v\equiv|{\bf v}|$ and the coefficients $X_{L}(v),X_{T}(v)$ and $Y(v)$ are equal to

[TABLE]

The equilibrium distribution function of the form

[TABLE]

with $T$ being the temperature of the plasma of light quarks and gluons, where heavy quarks are embedded, is expected to solve the transport equation (17). This is indeed the case if the coefficients $X^{ij}({\bf v})$ and $Y^{i}({\bf v})$ obey the condition

[TABLE]

which in the isotropic plasma reads

[TABLE]

When the plasma is isotropic and the coefficients $X_{L}(v)$ and $X_{T}(v)$ are equal to each other and independent of $v$, the Fokker-Planck equation reads

[TABLE]

where $X\equiv X_{L}(v)=X_{T}(v)$.

The quantities $X^{ij}({\bf v})$ and $Y^{i}({\bf v})$ have a clear physical meaning. As discussed in e.g. the classical monograph Kampen:1987 , the average momentum change per unit time and the correlation of momentum changes per unit time are given as

[TABLE]

Using the formulas (28) and (29), $X^{ij}({\bf v})$ and $Y^{i}({\bf v})$ can be related to the collisional energy loss $\frac{dE}{dx}$ and transverse momentum broadening $\hat{q}$ of a heavy-quark in the quark-gluon plasma, which play an important role in a theoretical description of the jet quenching phenomenon. The parameter $\hat{q}$ controls the radiative energy loss in a plasma medium Baier:1996sk . One easily finds that

[TABLE]

Since $\Delta x=v\Delta t$, the energy loss per unit path equals

[TABLE]

which in isotropic plasmas reads

[TABLE]

The coefficient $\hat{q}$, which is the broadening per unit path of the distribution of the test parton’s momentum transverse to the initial parton’s momentum, is immediately found as

[TABLE]

and in an isotropic plasma it equals

[TABLE]

When we deal with an equilibrium plasma and the coefficients $X_{L}(v),X_{T}(v)$ are equal to each other and independent of $v$, the Fokker-Planck equation can be related to the nonrelativistic Langevin equation Kampen:1987 . Then, the diffusion constant $D$ can be expressed as

[TABLE]

We are mostly interested in turbulent QCD plasmas populated with strong chromodynamic fields but we start with the equilibrium system where the fields are only at a level of thermal noise. We rederive the known Fokker-Planck equation Svetitsky:1987gq in a different way to demonstrate the reliability of our approach.

IV Equilibrium plasma

We assume that the quark-gluon plasma, in which heavy quarks are embedded, is in thermodynamical equilibrium and we first derive in this section explicit expressions for the coefficients $X_{L}$, $X_{T}$ and $Y$ which enter the Fokker-Planck equation.

IV.1 Computation of $X$ and $Y$

As the formula (16) shows, the quantity $X$ is given by the correlations functions $\langle E^{i}(t,{\bf r})E^{j}(t^{\prime},{\bf r}^{\prime})\rangle$, $\langle B^{i}(t,{\bf r})B^{j}(t^{\prime},{\bf r}^{\prime})\rangle$, $\langle E^{i}(t,{\bf r})B^{j}(t^{\prime},{\bf r}^{\prime})\rangle$, and $\langle B^{i}(t,{\bf r})E^{j}(t^{\prime},{\bf r}^{\prime})\rangle$ which were studied in detail in Mrowczynski:2008ae . The explicit expressions are collected in Appendix A. Since the correlation functions are of the structure (109), the tensor $X^{ij}$ is written as

[TABLE]

where the chromodynamic fields are expressed in the adjoint representation of the ${\rm SU}(N_{c})$ group and $\langle F^{i}_{a}F^{j}_{a}\rangle_{\omega,\,{\bf k}}$ is the fluctuation spectrum. For a translationally invariant system it is defined as

[TABLE]

A more general definition is discussed in Appendix B. Combining the equilibrium fluctuation spectra (110), (111) and (112), one finds

[TABLE]

where $\beta\equiv T^{-1}$ and $\varepsilon_{L,T}(\omega,{\bf k})$ are chromodielectric functions which for the equilibrium plasma of massless particles are also given in Appendix A.

After performing the elementary time integration in Eq. (36), one is left with the integral over the four-vector $k^{\mu}=(\omega,{\bf k})$. Taking into account only the terms of the integrand which are even as a function of $k^{\mu}$ and give nonzero contributions, one obtains

[TABLE]

where $\bar{\omega}\equiv{\bf k}\cdot{\bf v}$. In the limit $t\rightarrow\infty$, we have

[TABLE]

and thus we get

[TABLE]

One can show that the expression (41) properly approximates the formula (39) if the spectrum $\langle F^{i}_{a}F^{j}_{a}\rangle_{\omega,\,{\bf k}}$ weakly changes as a function of $\omega$ in the interval $[\bar{\omega}-\pi/t,\bar{\omega}+\pi/t]$. When the time grows the condition is easier and easier to fulfill.

Since the plasma under consideration is isotropic, the tensor $X^{ij}({\bf v})$ is fully determined by the two functions $X_{L}(v)$ and $X_{T}(v)$ introduced in Eq. (21). Substituting the explicit form of the fluctuation spectrum (38) into Eq. (41), one finds the coefficients $X_{L}(v)$ and $X_{T}(v)$ as

[TABLE]

where $C_{F}\equiv\frac{N_{c}^{2}-1}{2N_{c}}$ is the Casimir invariant.

The coefficient $Y^{i}({\bf v})$, which is determined by the correlations functions $\langle E^{i}(t,{\bf r})\delta Q_{0}({\bf r}^{\prime},{\bf p}^{\prime})\rangle$, and $\langle B^{i}(t,{\bf r})\delta Q_{0}({\bf r}^{\prime},{\bf p}^{\prime})\rangle$ can be derived directly from the formula (15). Such a derivation for a simplified case of only longitudinal electric field present in the plasma can be found in Mrowczynski:2009gf where it is shown that the condition (25) or (26) is indeed satisfied. Here instead we refer to the condition (25) to obtain $Y^{i}({\bf v})$.

Once the coefficients $X^{ij}({\bf v})$ is given by Eq. (21) together with Eqs. (42), (IV.1) and $Y^{i}({\bf v})$ by Eqs. (22) and (26), the Fokker-Planck equation (17) is fully specified. We note that $X^{ij}({\bf v})$ and $Y^{i}({\bf v})$ depend on heavy quark momentum ${\bf p}$ and its mass $m$ only through the velocity ${\bf v}={\bf p}/\sqrt{m^{2}+{\bf p}^{2}}$. Therefore, $X_{L}(v),~{}X_{T}(v)$ and $Y(v)$ become independent of $p$ when the heavy quarks of interest are truly relativistic and ${\bf p}^{2}\gg m^{2}$. Although, the coefficients $X_{L}(v),~{}X_{T}(v)$ and $Y(v)$ are independent of the quark mass, the Fokker-Planck equation (17) does depend on $m$ which is evident when the momentum derivatives are replaced by the velocity derivatives.

Since the coefficients $X_{L}(v),~{}X_{T}(v)$ and $Y(v)$ depend on the quark mass only through the velocity, one might think that the corresponding Fokker-Planck equation is applicable to quarks of any mass. However, it is not true. As mentioned in the introduction, the typical momentum transfer in a single collision, which is of order $gT$, must be much smaller than the quark momentum that is $gT\ll mv/\sqrt{1-v^{2}}$. And here the quark mass matters.

IV.2 Limit of small $v$

In the limit of small velocities of heavy quarks, the equilibrium Fokker-Planck equation gets a simpler form and the coefficients $X_{L}(v),~{}X_{T}(v)$ can be estimated analytically. Indeed, when heavy quarks are nonrelativistic (${\bf v}^{2}\ll 1$), we have $\bar{\omega}^{2}\ll{\bf k}^{2}$ and one can use the approximate formulas (116) and (117). Then, Eqs. (42) and (IV.1) read

[TABLE]

where the contributions originating from $\varepsilon_{T}(\omega,{\bf k})$ appear to vanish. Introducing spherical coordinates with the axis $z$ along the vector ${\bf v}$, the integrals (44) and (45) are rewritten as

[TABLE]

where the trivial azimuthal integrals are performed and $x\equiv\cos\theta$ with $\theta$ being the angle between ${\bf v}$ and ${\bf k}$.

When $T\gg m_{D}$, the integrals (46) and (47) can be estimated as follows. One first observes that the dominant contribution comes from $k\in[m_{D},T]$. Assuming that $\beta kxv\ll 1$, the integrals over $x$ are easily computed and one obtains

[TABLE]

Approximating the integrand as $k^{-1}$, we finally get

[TABLE]

As one sees in Eq. (48) or (49), $X_{L}(v)$ and $X_{T}(v)$ are independent of $v$ and equal to each other. The formula (35) is therefore applicable and the inverse diffusion constant equals

[TABLE]

which agrees with Eq. (12) of the study Moore:2004tg when only the logarithmic term is taken into account.

IV.3 Numerical results

Here we show some numerical results for the equilibrium QGP of $N_{c}=3$ and $N_{f}=2$. The Debye mass is computed according to the formula (115). Fig. 2 presents the coefficients $X_{L}(v)$ and $X_{T}(v)$, which are obtained directly from Eqs. (42) and (IV.1), as functions of the velocity $v$. The coupling constant and the temperature are $\alpha_{s}\equiv g^{2}/4\pi=0.1$ and $T=0.5~{}{\rm GeV}$. In Fig. 2 we show how $X_{L}(v)$ and $X_{T}(v)$ depend on the temperature $T$. The coupling constant is again $\alpha_{s}=0.1$ and the velocity equals $v=0.8$.

We have checked that the values of $X_{L}(v)$ and $X_{T}(v)$, which we obtained numerically, agree rather well with those computed by Svetitsky Svetitsky:1987gq except in the domain of small velocities $v\leq 0.1$. The agreement is not trivial because the coefficients of the Fokker-Planck equation were derived in Svetitsky:1987gq from the matrix elements of heavy quark binary interactions with plasma constituents. To remove infrared divergences of the matrix elements, a mass parameter, corresponding to the Debye mass, was included in the gluon propagator. Since the procedure is not very accurate, it presumably explains the difference with our $X_{L}(v)$ and $X_{T}(v)$ in the domain of small velocities. On the other hand our approach does not treat properly the interactions with a momentum transfer exceeding the Debye mass. Nevertheless the results of both approaches are numerically rather similar.

V Turbulent QGP

In this section we consider a Fokker-Planck equation of heavy quarks in a turbulent QGP which is populated with strong chromodynamic fields. The plasma is assumed to be isotropic and translationally invariant both in time and space. The tensor $X^{ij}({\bf v})$, which enters the Fokker-Planck equation (17), is given by Eq. (41). The method to obtain $Y^{j}({\bf v})$, which is used in Mrowczynski:2009gf , works only for equilibrium plasmas. Therefore, we will refer to the relation (25). However, it implicitly assumes that in the long time limit the system of heavy quarks described by the Fokker-Planck equation reaches a state of thermal equilibrium with temperature $T$. First of all, the value of $T$ is, in principle, unknown and one needs additional arguments to determine it. There is also a more important problem - it is unclear what are the properties of $X^{ij}({\bf v})$ for which the assumption of equilibrium makes sense. If, for example, there are only magnetic fields in the plasma, $X^{ij}({\bf v})$ is purely transverse and $TY^{i}({\bf v})=0$. Consequently, the Fokker-Planck equation reads

[TABLE]

and any stationary isotropic function $n(p)$ solves the equation because $\nabla_{p}n(p)\sim{\bf p}$. Therefore, pure magnetic fields do not drive the system to the thermal equilibrium, as expected. In spite of these concerns, we will use the relation (25) to get $Y^{j}({\bf v})$.

V.1 Gaussian correlation functions of independent $E$ and $B$ fields

We start with a simple model proposed in Asakawa:2006jn where the correlation functions of electric and magnetic fields are chosen in the following Gaussian form

[TABLE]

The correlation lengths $\sigma_{t},~{}\sigma_{r}$ and the parameters $M_{E},~{}M_{B}$ of dimension mass to the fourth power will be discussed later on. The fluctuation spectrum obtained from the correlation function (52) is also Gaussian and it equals

[TABLE]

Substituting the correlation functions (52), (53), and (54) into Eq. (41), the tensor $X^{ij}({\bf v})$ is found as

[TABLE]

which gives

[TABLE]

When $v^{2}\ll\sigma^{2}_{r}/\sigma^{2}_{t}$ and $v^{2}\ll M_{E}/M_{B}$, the coefficients $X_{L}(v)$ and $X_{T}(v)$ are, as in the equilibrium case, equal to each other and independent of $v$, and

[TABLE]

The Fokker-Planck equation is then of the form (27).

V.2 Gaussian correlation function of vector potentials

Since the electric and magnetic fields are, in general, coupled to each other, the functions $\langle E_{a}^{i}\,E^{j}_{b}\rangle$, $\langle E_{a}^{i}\,B^{j}_{b}\rangle$, $\langle B_{a}^{i}\,E^{j}_{b}\rangle$, and $\langle B_{a}^{i}\,B^{j}_{b}\rangle$ are not fully independent from each other. The electric and magnetic fields are automatically related to each other if one postulates the correlation function of the four-potential and then computes the correlation functions of the $E-$ and $B-$fields. In this section we follow this path. Specifically, we assume the Gaussian correlation function of the vector potential

[TABLE]

The parameter $M_{A}$ of the dimension mass squared will be discussed later on. The fluctuation spectrum of the potential equals

[TABLE]

We further choose the radiation gauge

[TABLE]

and the electric and magnetic fields are obtained in the linear regime as

[TABLE]

To get the fluctuation spectra $\langle E^{i}_{a}E^{j}_{b}\rangle_{\omega,{\bf k}}$, $\langle B^{i}_{a}E^{j}_{b}\rangle_{\omega,{\bf k}}$, $\langle B^{i}_{a}B^{j}_{b}\rangle_{\omega,{\bf k}}$ from the spectrum (61) we refer to the relation (126) derived in Appendix B. Using Eq. (63), the fluctuation spectra are found as

[TABLE]

Substituting the formulas (64), (65), and (66) into Eq. (41), one gets after some manipulations

[TABLE]

which gives

[TABLE]

When $v^{2}\ll\sigma^{2}_{r}/\sigma^{2}_{t}$, we find again, as in the equilibrium case, that the coefficients $X_{L}(v)$, $X_{T}(v)$ are equal to each other, but they do depend on $v$. Specifically,

[TABLE]

In contrast to the equilibrium result and Eq. (59), $X_{L}(v)$ and $X_{T}(v)$ vanish as $v\to 0$.

V.3 Stationary power spectrum of vector potential

When the momentum distribution of plasma constituents is anisotropic, the system is unstable due to the Weibel instability, see the review Mrowczynski:2016etf , and a strong chromomagnetic field is generated spontaneously. As shown in the numerical study Arnold:2005ef , the fluctuation spectrum of the soft fields becomes up to the wave number $k_{\rm max}$ stationary after a sufficiently long time and the spectrum decays with wave vector as $k^{-2}$. Inspired by these findings we choose, as previously, the radiation gauge (62) and consider the following spectrum

[TABLE]

where the parameters ${\cal M}$, $\mu$ and $k_{\rm max}$, which are all of the dimension of mass, will be determined later on. We note that $k_{\rm max}$ is of order of the Debye mass and that the small but nonzero parameter $\mu$ is introduced to eliminate the infrared divergence of the expression (71). In contrast to the study Arnold:2005ef , there are only zero frequency modes in the spectrum (71), and consequently the correlators involving electric field vanish. The only non-vanishing correlator is

[TABLE]

Using the formula (36), the tensor $X^{ij}({\bf v})$ is found as

[TABLE]

The coefficient $X_{L}(v)$ given by Eq. (23) vanishes because of transversality of magnetic field while $X_{T}(v)$ equals

[TABLE]

After taking the elementary integrals, one obtains

[TABLE]

As seen, $X_{T}(v)$ vanishes when $v\to 0$ and the formula (75) simplifies to

[TABLE]

when $k_{\rm max}\gg\mu$.

In Tables 1 and 2 there are collected the coefficients $X_{L}(v)$, $X_{T}(v)$ and corresponding transport coefficients $\frac{dE}{dx}$, $\hat{q}$ obtained in the models of turbulent plasma discussed in the subsections V.1, V.2 and V.3. The models are called ‘Gauss E & B’, ‘Gauss A’ and ‘Stationary A’, respectively. We observe that $\frac{dE}{dx}$ is not always reliable because there is no guarantee that heavy quarks evolve in turbulent plasmas towards the state of thermodynamical equilibrium as required by the relation (26) that is used. We also note that the formula of $\hat{q}$ in the ‘Gauss E & B’ model, which holds not for a test quark with any $v$ but for a test gluon with $v=1$, was first obtained in Majumder:2007zh .

VI Equilibrium vs. isotropic turbulent plasma

We are going to compare the equilibrium plasma to the turbulent one at the same energy density. The energy density of the equilibrium quark-gluon plasma of $N_{f}$ massless flavors equals

[TABLE]

The density of energy accumulated in chromodynamic fields is expressed as

[TABLE]

If the fluctuation spectra $\langle E^{i}_{a}E^{i}_{a}\rangle_{\omega,\,{\bf k}}$ and $\langle B^{i}_{a}B^{i}_{a}\rangle_{\omega,\,{\bf k}}$ are of the Gaussian form (55), the energy density equals

[TABLE]

When the fluctuation spectra $\langle E^{i}_{a}E^{i}_{a}\rangle_{\omega,\,{\bf k}}$ and $\langle B^{i}_{a}B^{i}_{a}\rangle_{\omega,\,{\bf k}}$ are given by Eqs. (64) and (65), we have

[TABLE]

Finally, if the fluctuation spectrum $\langle A^{i}_{a}A^{i}_{a}\rangle_{\omega,\,{\bf k}}$ is of the form (71), the purely magnetic energy density equals

[TABLE]

where we assume that $k_{\rm max}\gg\mu$. We note that in the weak coupling limit the magnitudes of electric and magnetic fields in turbulent plasmas, which are according to the above estimates of the order $T^{2}$, are much larger than in equilibrium plasmas where the fields are typically of the order $g^{2}T^{2}$. This is the main reason why the turbulent plasmas, which are discussed here, are qualitatively different than the equilibrium one.

To reduce the number of parameters of our models of turbulent plasma we adopt the simplifying assumptions that $M_{E}=M_{B}$ and $\sigma_{t}=\sigma_{r}$. Then, when the energy density of the fields in turbulent plasma (79), (80) and (81) equals the energy density of the equilibrium QGP (77), the parameters $M_{E}$, $M_{A}$ and ${\cal M}$ are equal to

[TABLE]

We set $N_{c}=3$ and $N_{f}=2$. With the simplifying relations $M_{E}=M_{B}$ and $\sigma_{t}=\sigma_{r}$ and the equalities (82), one obtains the parametric estimates of the coefficient $X_{T}(v)$ which are shown in the second column of Table 3 where we have additionally assumed that $v^{2}\ll 1$.

Since the parameters $\sigma_{t},~{}\sigma_{r}$ determine the field correlation length in time and space, they are of order of the screening length and thus we choose $\sigma_{t}=\sigma_{r}=m_{D}^{-1}$. Similarly, $k_{\rm max}$ is a soft momentum and following Arnold:2005ef we set $k_{\rm max}=5m_{D}$. Using the formula (115), we get

[TABLE]

The parametric estimates of $X_{T}(v)$, which take into account the relations (83), are shown in the third column of Table 3. As seen, the interaction of heavy quarks is much stronger in the turbulent plasma than in the equilibrium one if the plasma is truly weakly coupled with $g\ll 1$. The effect disappears when the coupling constant $\alpha_{s}\equiv g^{2}/4\pi$ is of realistic value and gives $g$ close to unity or even bigger.

One asks whether the strong magnetic fields discussed above are stable or may be they rapidly decay due to the Nielsen-Olesen instability Nielsen:1978rm which has been discussed in the context of quark-gluon plasma in Berges:2011sb ; Kurkela:2012hp . Strictly speaking, the Nielsen-Olsen instability occurs in systems with a homogeneous magnetic field because the system’s energy can be reduced by particles circulating in the magnetic fields. The unstable mode grows as $e^{\gamma t}$ with $\gamma=\sqrt{gB}$. When the field is inhomogeneous with the wave vector $k$, the decrement of growth is reduced to $\gamma=\sqrt{gB-k^{2}}$. The instability disappears when $gB<k^{2}$ which physically means that the length of inhomogeneity is shorter than the Larmor radius. As discussed above, $B\sim T^{2}$ and $k\sim gT$ in the turbulent plasma and thus the Nielsen-Olsen instability is unavoidable as long as $g\ll 1$. The lifetime of strong magnetic is then of order $(\sqrt{g}\,T)^{-1}$. In Sec. VII we discuss the glasma Gelis:2012ri which is of our main interest. We show that the instability is absent in the glasma because the field correlation length is of the order of the inverse stauration scale $Q_{s}$.

In Figs. 4 and 4 we present how $X_{L}(v)$ and $X_{T}(v)$ depend on the velocity $v$ in the plasma models under consideration. The coupling constant and temperature are $\alpha_{s}=0.1$ and $T=0.5~{}{\rm GeV}$. Figs. 6 and 6 show the corresponding energy loss and momentum broadening as functions of $v$. Since the energy density of the turbulent plasma is fixed, the temperature $T$, which enters Eq.  (25) or (26), is fixed as well. The coefficient $\hat{q}$ is computed for $v\geq 0.3$ because it diverges as $v\to 0$ in the model ‘Gauss E & B’. However, $\hat{q}$ is of no physical relevance for small velocities. As seen, the coefficients $X_{L}(v)$ and $X_{T}(v)$ and consequently $dE/dx$ and $\hat{q}$ are rather different in the equilibrium and turbulent plasmas - not only the magnitudes differ but the dependences on the heavy quark’s velocity are different. The interaction of heavy quarks is particularly strong in the model ‘Gauss E & B’.

One wonders why the models ‘Gauss E & B’ and ‘Gauss A’, which at first glance are quite similar, give rather different results shown Figs. 4 and 4. Comparing the fluctuation spectra (64) and (65) to (55), one observes that the low frequency and long wavelength fields are suppressed in the ‘Gauss A’ model. These fields appear to contribute more effectively to the transport coefficients than the high frequency and short wavelength fields.

VII Glasma

When relativistic heavy ions collide color charges of partons confined in the colliding nuclei generate strong chromodynamic fields right after the collision. Since the system of infinitely contracted nuclei moving against each other with the speed of light is boost invariant, so is the configuration of generated fields. As shown in the detailed analytic study Chen:2015wia , the chromoelectric and chromomagnetic fields spanned between the receding nuclei are initially only parallel to the beam direction. Transverse field components start to develop later on. We focus here on the longitudinal fields which dominate the glasma’s dynamics and are invariant under Lorentz transformations along the beam direction identified with the axis $z$.

VII.1 Field correlation functions

Since the electric and magnetic fields are expressed in the Abelian limit through the four-potential according to Eq. (63), the potential generating the fields only along the axis $z$ is of the form

[TABLE]

that is the [math] and $z$ components of $A^{\mu}$ depend on $t$ and $z$ while the $x$ and $y$ components on $x$ and $y$. The electric field corresponding to the potential (84) depends on $t$ and $z$ while the magnetic field on $x$ and $y$. The electric field $E^{z}_{a}$ is boost invariant if it depends on $t$ and $z$ only through the proper time $\tau\equiv\sqrt{t^{2}-z^{2}}$ which is the Lorentz scalar. However, we do not require the boost invariance of the fields but of the field correlators. To write down a general expression of the electric field correlator $\langle E^{z}_{a}(t_{1},z_{1})\,E^{z}_{b}(t_{2},z_{2})\rangle$, we introduce the variables

[TABLE]

and we note that the proper times $\tau_{i}$ and space-time rapidities $\eta_{i}$ are well defined only for the time-like two-vectors $(t_{i},z_{i})$. The boost invariant correlation function of the electric fields is assumed to be

[TABLE]

where $\Theta(t_{i}^{2}-z_{i}^{2})$ is the Heaviside step function and $f_{E}(\tau,\eta)$ is an arbitrary function. Since $\tau_{i}$ and $\eta_{i}$ are well defined only for $t_{i}^{2}\geq z_{i}^{2}$, we require that the correlation function (86) vanish when the space-time points $(t_{i},z_{i})$ are localized beyond the light cone of the point $(0,0)$. We observe that the correlation function, which depends not only on $\tau_{1}-\tau_{2}$ but on both $\tau_{1}$ and $\tau_{2}$ , is also boost invariant but we assume the translation invariance in the $\tau$ variable just for simplicity.

The boost invariant magnetic field correlator is chosen to be

[TABLE]

that is the plasma is assumed to be translationally invariant in the $x$ and $y$ directions. We note that the fields $E^{z}_{a}(t,z)$ and $B^{z}_{a}(x,y)$ are completely decoupled from each other in the Abelian limit and thus the mixed correlator $\langle E^{z}_{a}(t_{1},z_{1})\,B^{z}_{b}(t_{2},z_{2})\rangle$ is expected to be small or vanish.

We introduce the unit vector ${\bf n}=(0,0,1)$ along the axis $z$ and the correlation functions of electric and magnetic fields are written as

[TABLE]

The correlators are chosen to be of the Gaussian form

[TABLE]

with the real positive parameters $\tilde{M}_{E}$, $\tilde{M}_{B}$, $\sigma_{\tau}$, $\sigma_{\eta}$ and $\sigma_{T}$ to be determined later on.

VII.2 Computation of $X$ and $Y$

Substituting the correlation function (89) with (91) into Eq. (16), the magnetic contribution to the tensor $X^{ij}({\bf v})$ is found as

[TABLE]

where

[TABLE]

Since the Fokker-Planck equation is derived in the long time limit, we assume that $t\gg\sigma_{T}$ and the integral becomes Gaussian. Thus, we get

[TABLE]

where $v_{T}\equiv|{\bf v}_{T}|=\sqrt{v_{x}^{2}+v_{y}^{2}}$ is the quark velocity perpendicular to ${\bf n}$.

The electric contribution to the tensor $X^{ij}({\bf v})$ is

[TABLE]

where

[TABLE]

The integral from Eq. (95) with the correlation function (90) is difficult to compute analytically. The problem greatly simplifies in the long time limit. When $t\gg z$ and $t\gg\sigma_{\tau}$ one easily finds

[TABLE]

which is independent of $\sigma_{\eta}$. We have checked numerically a quality of the approximation (98). In Fig. 7 we show the numerically computed coefficient $X^{zz}_{E}({\bf v})$ divided by the approximate expression (98) for $t=10\sigma_{\tau}$ and $t=30\sigma_{\tau}$ as a function of heavy-quark velocity $v_{z}$. The coefficient is independent of $v_{x}$ and $v_{y}$. The correlation length in rapidity is chosen as $\sigma_{\eta}=1$ but for a bigger $\sigma_{\eta}$ the approximation (98) is even more accurate. We conclude that the approximation (98) works pretty well and we write the tensor $X^{ij}({\bf v})$, which includes the magnetic and electric contributions, as

[TABLE]

Using the relation (25), one finds that the vector $Y^{i}({\bf v})$ equals

[TABLE]

The temperature will be estimated in the next section but, as already discussed at the beginning of Sec. V, the result (100) should be treated with a reservation.

VII.3 Estimates of parameters

To get approximate values of the parameters, which characterize the glasma, let us estimate the density of energy released in relativistic heavy-ion collisions. When one deals with the central collisions of two nuclei of mass number $A$, the energy density in the center-of-mass frame of colliding nuclei is roughly

[TABLE]

where $\sqrt{s}$ is the center-of-mass energy of nucleon-nucleon collision, $c_{\rm inel}$ is the inelasticity parameter which gives the fraction of all accessible energy going to particle production, $R_{A}$ is the radius of colliding nuclei and $l$ is a length of the cylinder where the energy is released. Assuming that $c_{\rm inel}=0.5$ Navarra:2003am and taking $A=200$, $R_{A}=7$ fm and $l=1$ fm, one obtains $\varepsilon_{\rm coll}\approx 3.25\;{\rm TeV/fm^{3}}$ for $\sqrt{s}$ = 5 TeV which is the energy of Pb-Pb collisions at LHC of the 2015 run.

According to Eq. (78) the density of energy accumulated in chromodynamic fields equals

[TABLE]

where the explicit form of the correlation functions (90) and (91) have been used. The field energy density is controlled by the parameters $\tilde{M}_{E}$ and $\tilde{M}_{B}$ but is independent of the correlations lengths $\sigma_{T}$, $\sigma_{\tau}$ and $\sigma_{\eta}$. Assuming $\tilde{M}_{E}=\tilde{M}_{B}$, as suggested Chen:2015wia , we get $\varepsilon_{\rm field}=8\tilde{M}_{E}$ for $N_{c}=3$. Requiring that $\varepsilon_{\rm field}=\varepsilon_{\rm coll}$, the parameter $\tilde{M}_{E}$ is estimated as

[TABLE]

To use the formula (25) to obtain the coefficient $Y^{i}({\bf v})$, one needs a temperature of the equilibrated quark-gluon plasma of the same energy density as the glasma. Using the formula (77), the equation $\varepsilon_{\rm field}=\varepsilon_{\rm QGP}$ provides

[TABLE]

where $N_{c}=3$ and $N_{f}=2$.

Within the CGC approach the strong longitudinal chromodynamic fields are screened on transverse distance which is of the order of the inverse saturation momentum $Q_{s}$. Since $Q_{s}$ is estimated as 2 GeV Gelis:2012ri , we choose the transverse correlation $\sigma_{T}=Q_{s}^{-1}=0.5~{}{\rm GeV}^{-1}$. We also assume that $\sigma_{\tau}=\sigma_{T}$. The remaining parameter is the coupling constant which, as previously, is chosen to be $\alpha_{s}=0.1$. As already mention in Sec. VI, the glasma magnetic field is stable against the Nielsen-Olesen instability, because the inhomogeneity length of order $Q_{s}^{-1}$ is not much bigger than the Larmor radius $(gB)^{-1/2}$, $Q_{s}$ is even bigger than $\tilde{M}_{B}^{1/4}$.

Using the formulas (100) and (99) combined with Eqs. (31) and (34), we obtain the following estimates of the energy loss and momentum broadening

[TABLE]

where $\theta$ is the angle between ${\bf v}$ and ${\bf n}$. Because the electric field is along the beam axis, the collisional energy loss is maximal when the heavy quark moves along the axis, ${\bf v}\parallel{\bf n}$. The maximal momentum broadening occurs when ${\bf v}\perp{\bf n}$.

Let us compare the numerical values (105) and (106) with those which are required to properly model experimental data on the charm meson suppression. The collisional energy loss and momentum broadening of a charm quark with 10 GeV momentum in the plasma of the temperature from the interval 0.35 – 0.5 GeV are estimated Prino:2016cni as

[TABLE]

As seen, the values (105) and (106) can be significantly larger than (107) and (108), suggesting that in spite of a short lifetime of the glasma it can provide a significant contribution to the collisional and radiative energy loss to heavy quarks from relativistic heavy-ion collisions. Consequently, the effect should be included in the phenomenology of jet quenching.

VIII Summary, Conclusions and Outlook

Applying the so-called quasi-linear theory Ved61 ; LP81 , we have derived a general form of the Fokker-Planck equation of heavy quarks embedded in the plasma of light quarks and gluons. Since the interaction is taken into account through the correlation functions of chromodynamic fields, heavy quarks are seen as interacting not with plasma constituents but rather with fields present in the plasma. At first we have obtained the explicit form of the equation for the case of equilibrium plasma which was studied long ago Svetitsky:1987gq using the standard method where the Fokker-Planck equation simply approximates the Boltzmann one. Although our approach is noticeably different, the Fokker-Planck equation we obtained agrees with the standard one Svetitsky:1987gq .

In the second part of the paper, the method developed for the equilibrium plasma has been applied to the turbulent plasma populated with strong fields. The parametric estimate shows that the interaction of heavy quarks with the turbulent plasma is much stronger than with the equilibrium one of the same energy density if the coupling constant is truly small. The effect is less prominent for a realistic value of the coupling constant and the difference depends on characteristics of the plasma fields. Within the ‘Gaussian E & B’ model both the energy loss and momentum broadening are significantly bigger than the equilibrium results. A dependence of $dE/dx$ and $\hat{q}$ on heavy quark velocity also strongly depends on how the turbulent plasma is modeled.

The third part of our study is devoted to the glasma from the earliest stage of relativistic heavy-ion collisions. Assuming that there are chromoelectric and chromomagnetic fields only along the beam direction we have derived the appropriate Fokker-Planc equation. We have also shown that in spite of its short lifetime the glasma can provide a significant contribution to the collisional and radiative energy loss of heavy quarks.

Our findings clearly suggest a direction of further work. We need a more realistic model of turbulent QCD plasma from relativistic heavy-ion collisions. In contrast to the simple model discussed here, a temporal evolution of the glasma has to be taken into account and the fields cannot be purely longitudinal. The CGC studies Gelis:2012ri and, in particular, the analytic analysis Chen:2015wia provide a very good guidance to build up such a model.

Acknowledgments

I am very grateful to Margaret Carrington for numerous fruitful discussions and valuable comments.

Appendix A Equilibrium correlation functions

The correlation functions of chromodynamic fields in the equilibrium plasma, which were studied in detail in Mrowczynski:2008ae , can be expressed as

[TABLE]

where $H^{i}_{a}(t,{\bf r})$ and $K^{j}_{b}(t,{\bf r})$ is either the electric or magnetic field, and the fluctuation spectra are

[TABLE]

$\epsilon^{ijm}$ is the antisymmetric tensor, $\beta\equiv T^{-1}$ with $T$ being the system’s temperature and $\varepsilon_{L,T}(\omega,{\bf k})$ are chromodielectric functions. For an equilibrium plasma of massless particles, the functions are well-known to be, see e.g. lebellac ,

[TABLE]

with $m_{D}$ being the Debye mass which for the quark-gluon plasma of $N_{f}$ massless flavors equals

[TABLE]

When $\omega^{2}\ll{\bf k}^{2}$, the dielectric functions can be approximated as

[TABLE]

We note that the fluctuation spectra of pure classical fields were actually derived in our study Mrowczynski:2008ae . The effect of Bose statistics of field quanta has been included in the formulas (110), (111), and (112) by means of the substitution

[TABLE]

The absolute value of the frequency results from the following reasoning. Since the electric and magnetic fields are real in coordinate space, their correlation functions are real as well. Consequently, the fluctuation spectra must obey

[TABLE]

One checks that the formulas (110), (111), and (112) indeed satisfy the symmetry (119) and thus the correlation function (109) is real as it should be.

Appendix B Correlation function and fluctuation spectrum

We discuss here the relation between the correlation function of the Fourier transformed fields and the fluctuation spectrum. The correlation function of the Fourier transformed fields equals

[TABLE]

where $H^{i}_{a}(t,{\bf r})$ and $K^{j}_{b}(t,{\bf r})$ is either the electric field, magnetic field, or the potential. To define the fluctuation spectrum we first write down the correlation function as

[TABLE]

where the new space-time variables read

[TABLE]

The fluctuation spectrum is defined as

[TABLE]

In general, the fluctuation spectrum (124) depends on $t$ and ${\bf r}$. However, if the system is stationary and homogeneous, that is translationally invariant in both space and time, the spectrum is independent of $t$ and ${\bf r}$. Therefore, Eq. (124) can be rewritten as

[TABLE]

where $V{\cal T}$ is the space-time volume occupied by the system. Substituting the space-time correlation function expressed through $\langle H^{i}_{a}(\omega,{\bf k})\,K^{j}_{b}(\omega^{\prime},{\bf k}^{\prime})\rangle$ into Eq. (125) and performing the trivial integrations involving delta functions, one obtains the desired relation

[TABLE]