An Effective Field Theory Study of Medium Heavy Quark Evolution

Abstract

The evolution of hard probes in a medium is a complex multiscale problem that significantly benefits from the use of Effective Field Theories (EFTs). Within the EFT framework, we aim to define a series of EFTs in a way that addresses each energy scale individually in separate steps. However, studying hard probes in a medium presents challenges. This is because an EFT is typically constructed by formulating the most general Lagrangian compatible with the problem’s symmetries. Nevertheless, medium effects may not always be encoded adequately in an effective action. In this paper, we construct an EFT that is valid for studying the evolution of a heavy quark in a QCD plasma containing few other heavy quarks, where degrees of freedom with an energy of the order of the temperature scale are integrated out. Through this example, we explicitly demonstrate how to handle the doubling of degrees that arise in non-equilibrium field theory. As a result, we derive a Fokker–Planck equation using only symmetry and power counting arguments. The methods introduced in this paper will pave the way for future developments in the study of quarkonium suppression.

1. Introduction

Quark–gluon plasma is a state of matter that forms at high temperatures and densities, in which quarks and gluons are not confined within hadrons. This state of matter can be created on Earth in experiments using ultrarelativistic heavy-ion collisions. However, quark–gluon plasma exists only for a very short time during these collisions, so we must study the particles produced during this brief period to learn about the properties of the plasma. One approach is to study “hard probes”, which are observables that are both significantly affected by the medium and can be measured in the challenging environment of a heavy-ion collision. Examples of hard probes include heavy quarks, heavy quarkonium, and jets.

Studying heavy particles and jets in a medium requires dealing with largely separated energy scales. Particles with energies of the order of the heavy quark mass or a hard parton energy are rare in the medium and can be accurately described using perturbative QCD. However, particles with energies of the order of the temperature are sensitive to the medium and cannot always be described using perturbation theory. Therefore, it is interesting to separate the contribution from these two kind of particles to encode non-perturbative effects in parameters or functions that can be computed using non-perturbative methods like lattice QCD. It is also important to be careful when describing systems with largely separated energy scales. On one hand, the appearance of largely separated energy scales can lead to a breaking of naive perturbation theory, meaning that the size of a contribution cannot be directly related to its number of loops. On the other hand, non-perturbative studies using lattice QCD are also challenging because a very large lattice is required to accommodate all the energy scales.

These problems can be solved by using EFTs. An EFT is a quantum field theory that gives the same results as another more general theory at low energies. They are constructed in the following way [1]:

• Identify the relevant degrees of freedom and the symmetries of the problem.
• Write the more general Lagrangian that respects the symmetries of the problem using the relevant degrees of freedom.
• An EFT must be equipped with a power counting. This means that there is a simple rule to predict how large the contribution is of each term in the Lagrangian for a given observable.
• The unknown parameters in the effective Lagrangian are called Wilson coefficients. They are fixed by imposing that the EFT gives the same results as the full theory at low energies. This procedure is called matching.

Note that an EFT Lagrangian has an infinite number of terms. However, the theory still has predictive power thanks to the power counting.

The use of EFTs can also allow relating physical observables with quantities computable on the lattice in a more direct way. A good example is the description of heavy quarkonium. Heavy quarkonium is a bound state made of a heavy quark and a heavy antiquark. Quarkonium is a non-relativistic system in which well-separated energy scales appear. They are the mass, M, the inverse of the typical radius, $1/r\sim Mv$ with $v\ll 1$, and the binding energy $E\sim M{v}^2$. Non-relativistic QCD (NRQCD) [2,3] is an EFT valid for energies smaller than M. Since $M\gg {\mathsf{\Lambda }}_{QCD}$, we can reliably match QCD to NRQCD using perturbation theory. We can also use Potential NRQCD (pNRQCD) [4,5,6], an EFT valid for energies much smaller than $1/r$. Since this scale is not always perturbative, there are cases in which the matching between NRQCD and pNRQCD cannot be done in perturbation theory. However, since each term in the pNRQCD Lagrangian has a scaling with $1/M$ that is easy to determine, we can perform the matching between NRQCD and pNRQCD in the limit in which the mass of the heavy quarks is infinite (the static limit). Once the matching is carried out, we can use these Wilson coefficients to compute the properties of quarkonium states with finite heavy quark masses. This is a rigorous way of showing that the spectrum of quarkonium can be determined by solving a Schrödinger equation in which the potential is computed on the lattice using static quarks.

EFTs have been applied to the study of hard probes in heavy-ion collisions. In the case of heavy quarkonium, medium modified versions of NRQCD and pNRQCD have been used to compute the thermal corrections to the mass and the medium-induced decay width [7,8,9]. pNRQCD has also been used to study the evolution of the population of bound states inside of a medium in the cases in which $1/r\gg T$ [10,11,12,13]. Regarding jets, there are recently developed EFTs for the study of jet broadening in [14] as well as jet sub-structures in heavy-ion collisions [15].

However, the application of the EFT framework to the study of hard probes in a medium presents a conceptual challenge. The procedure that we have outlined before needs to be modified when the medium affects the matching between the full theory and the EFT. At a finite temperature, computations of non-static properties are performed using so-called real-time formalism [16]. In this formalism, there is a doubling of the degrees of freedom. This means that the standard path integral has to be substituted by a path integral in which the time integration follows the complex Schwinger-Keldysh contour (see Figure 1). In practice, instead of considering explicitly the complex–time integration, we name fields with a time argument on the upper (lower) branch of the contour fields of type 1(2). As an example, let us consider a QFT describing the evolution of a field $\varphi$. Then, the action can be written as

$$S={\int }_{C}dtL\left[\varphi \left(t\right)\right]={\int }_{-\infty }^{\infty }dtL\left[{\varphi }_1\left(t\right)\right]-{\int }_{-\infty }^{\infty }dtL\left[{\varphi }_2\left(t\right)\right],$$

(1)

where the subindex C means that the integration is performed along the Schwinger-Keldysh contour. In this way, due to the properties of the path integral, correlators in which only fields of type 1(2) appear are chronologically (anti-chronologically) ordered. However, we might be interested in evaluating correlators that are neither chronologically nor anti-chronologically ordered. For example, if we want to gain access to the distribution function, we would need to compute $\langle \varphi \left(t\right)\varphi \left(0\right)\rangle$, where the ordering of the operators is as written. Note also that, although in Equation (1) there are no terms containing, at the same time, both types of fields, the propagator of the field ${\varphi }_i$ is not diagonal in the $1,2$ indices [16].

What we have explained up until now regarding the Schwinger-Keldysh contour (for the specific case of a QFT describing the evolution of a field $\varphi$) is standard textbook material. However, it is important to recall this information in order to highlight non-trivial features that appear when we apply the EFT framework to a non-equilibrium or thermal field theory. When we integrate out degrees of freedom, the matching between the full theory and the EFT might induce terms including the two types of fields. In other words, the EFT might not show the structure that appears in Equation (1). Physically, we might interpret this in the following way. The EFT describes the evolution of an open quantum system interacting with an environment (the medium degrees of freedom that have been integrated out). It is a well-known result in quantum mechanics that the evolution of an open quantum system cannot be encoded in a Hamiltonian [17]. In other words, in general, an operator does not exist, such that the time-evolution of the density matrix can be written as a von-Neumann equation. Instead, the evolution of the density matrix of a open quantum system follows a more general kind of equation. The same is true if we apply path integral formalism. The effect of high-energy fields on low-energy ones can not be encoded into an effective action when we are dealing with an open quantum system. We need a more general object, called a influence functional, instead [18]. For the example of a $\varphi$ field

$$\langle \mathcal{P}{\varphi }_i\left(t\right)\varphi \left(t^{\prime }\right)\rangle =\int \mathcal{D}{\varphi }_1\mathcal{D}{\varphi }_2\mathcal{F}[{\varphi }_1,{\varphi }_2]{\varphi }_i\left(t\right){\varphi }_j\left(t^{\prime }\right),$$

(2)

where $\mathcal{P}$ means ordering along the Schwinger-Keldysh contour and $\mathcal{F}$ is the influence functional. Only in the case of an isolated system does $\mathcal{F}={\mathrm{e}}^{S\left[{\varphi }_1\right]-S\left[{\varphi }_2\right]}$, where S is the effective action.

Taking this into account is a substantial modification of the second point of the procedure we outlined to construct an EFT. Up until now, this difficulty has been skipped in the study of hard probes using two different strategies. The first one is to focus on computing the binding energy and the decay width of non-relativistic systems. In this case, the effects of the doubling of the degrees of freedom are neutralized in the approximation in which the heavy particles are diluted [7,8]. The second strategy is to study cases in which the temperature is smaller than the energies that are integrated out to define the EFT [10,11].

The construction of EFTs that do not consist in an effective action has been investigated in the context of the development of an EFT for hydrodynamics [19,20,21,22]. However, to our knowledge, these developments have not been applied to the study of hard probes in a QCD plasma. In this manuscript, we aim to fill this gap by studying a simple example: a heavy quark close to thermalization in a medium with few other heavy quarks. In this problem, we can identify three widely separated energy scales.

• The mass of the heavy quark M.
• The spatial momentum of the heavy quark p, which is of the order of $\sqrt{TM}$, with T being the temperature.
• The temperature T and other energy scales that the medium might induce.

It is well known that this situation can be studied using a Fokker–Planck equation or the physically equivalent Langevin equation [23,24,25,26,27,28,29]. In this manuscript, we will give an EFT perspective on the problem. Following the EFT philosophy, we will deal with each energy scale in a separate way. We will follow the following steps:

• As a starting point, we can use NRQCD. This automatically encodes the effects of the degrees of freedom with an energy of order M.
• It is straight-forward to define an EFT for energies smaller than $\sqrt{MT}$. This EFT is a modification of NRQCD using the momentum label technique common in Soft-Collinear Effective Theory (SCET) [30,31] and On-Shell Effective Theory (OSEFT) [32,33]. An equivalent EFT, with the same degrees of freedom, symmetries, and power counting, was previously introduced to study a completely different problem. This is an NRQCD for semi-hard fields (${\mathrm{NRQCD}}_{sh}$) [34]. This EFT was originally introduced to study modes with an energy of the order $\sqrt{M{\mathsf{\Lambda }}_{QCD}}$ in the study of heavy quarkonium at $T=0$. In our case, the role of ${\mathsf{\Lambda }}_{QCD}$ is played by T. However, symmetries and power counting arguments remain the same.
• We define an EFT for energies below T. In this EFT, there are terms that mix fields of type 1 and type 2. These are dealt with following the approach of [19] but supplementing it with the properties of the dilute expansion. We call this EFT Langevin Effective Theory (LET) due to its connection with the Langevin equation.
• We perform the matching between ${\mathrm{NRQCD}}_{sh}$ and LET in the one-gluon exchange approximation as an illustration. This allows us to encode the influence of the degrees of freedom with an energy of order T.
• We use LET to compute the evolution of the density of heavy quarks, obtaining a Fokker–Planck equation for the density that coincides with the evolution resulting from a Langevin equation.

This procedure could seem extremely complex to obtain results that are already well-known. However, our aim is to develop techniques that could be useful for more interesting cases. In particular, we have in mind the study of the evolution of the reduced density matrix of quarkonium in the regime $T\sim \frac{1}{r}$. The problem that we study in this manuscript is simpler than quarkonium in the regime $T\sim \frac{1}{r}$. However, in both cases we are dealing with dilute non-relativistic particles. One important observation that the framework developed here makes manifest is that the physics at the scale T can be studied in the static limit. In other words, the matching between ${\mathrm{NRQCD}}_{sh}$ and LET can be done in the static limit, and we can use those Wilson coefficients to study the finite M case. A similar thing happens regarding the matching between NRQCD and pNRQCD at $T=0$. These two results combined strongly support the idea that the matching between NRQCD and pNRQCD at a finite temperature can be also carried out in the static limit. We note that the evolution of a static heavy quark–antiquark pair in terms of gauge-invariant expectation values has already been discussed in [35].

As a complement to the previous study, we will also discuss the case in which $p\sim T$. In this case, we integrate out the scales p and T at the same time. In other words, we go directly from NRQCD to LET. We observe that the structure of LET in this case is still the same as in the case $p\gg T$. However, the Wilson coefficients are different. We observe that there is a smooth transition between the cases $p\sim T$ and $p\gg T$, as should be the case. Physically, if we start with a heavy quark at rest, it will start to gain momentum due to broadening. As the momentum becomes larger, two complementary things happen. On the one hand, the Fokker-Planck equation becomes accurate. On the other hand, the drag force becomes a leading order effect.

The present work has some similarities with [36]; however, the perspective is quite different. In [36], the AdS/CFT correspondence in Supersymmetric Yang–Mills is taken as a starting point, while, in this work, we tried to be as agnostic as possible about the properties of the medium.

The manuscript is organized as follows. In Section 2, we will discuss ${\mathrm{NRQCD}}_{sh}$. Next, in Section 3 , we will present LET. Section 4 discusses the computation of the evolution of the distribution of heavy quarks within LET. In Section 5, we discuss the matching between ${\mathrm{NRQCD}}_{sh}$ and LET. In Section 6, we consider the case of a heavy quark with momentum $p\sim T$. Finally, in Section 7, we give our conclusions.

2. ${\mathbf{NRQCD}}_{\mathit{sh}}$

Our starting point is the NRQCD Lagrangian [2,3]. To fix the notation, we write its heavy quark sector

$${\mathcal{L}}_{\psi }={\psi }^{†}\left(x\right)\{i{D}_0+\frac{{\mathbf{D}}^2}{2M}+c_4\frac{{\mathbf{D}}^4}{8M^3}+c_{F}g\frac{\sigma ·\mathbf{B}}{2M}\}\psi \left(x\right)+\cdots$$

(3)

Here, $D_{\mu }={\partial }_{\mu }-ig{A}_{\mu }$, $\sigma$ is a Pauli matrix and E and B are the chromoelectric and chromomagnetic fields. We use the following power counting. The spatial momentum of the heavy quark scales like $\sqrt{MT}$ plus a possible residual momentum of order T. The gauge field is only sensitive to the scale T. In Equation (3), we have only written terms of the order of ${\psi }^{†}\psi \frac{T^2}{M}$ or lower. These are the terms that we need in order to compute $\frac{1}{f}\frac{df}{dt}$ to order $\frac{T^2}{M}$ where $f(p,\mathbf{R})$ is the distribution function of heavy quarks. However, only the first two terms in the Lagrangian will end up contributing. We define f in terms of the < propagator of heavy quarks. More details will be given later.

Using NRQCD to study this problem is not completely optimal. As we mentioned before, the spatial momentum of heavy quarks is of the order $\sqrt{MT}$ while each interaction with the medium changes the momentum of the heavy quark by an amount of the order T. This implies that the term ${\psi }^{†}\frac{{\mathbf{D}}^2}{2M}\psi$ hides contributions of different sizes. Then, the power counting is not completely clear. Each time that we apply this term it is not obvious whether we will get a contribution of order T, $T\sqrt{\frac{T}{M}}$ or $\frac{T^2}{M}$. At best, we can put an upper bound on the size of the contribution. To improve the situation, we can introduce momentum-label fields, as is performed in SCET and OSEFT. Let us divide the spatial momentum of heavy quarks into two pieces

$$\mathbf{p}+\mathbf{k},$$

(4)

where $\mathbf{p}$ is of the order $\sqrt{MT}$ and $\mathbf{k}$ is a residual momentum of the order T. To improve this, we can perform the following transformation:

$$\psi (t,\mathbf{x})=\sum _{\mathbf{p}}{\mathrm{e}}^{i\mathbf{p}\mathbf{x}}{\xi }_{\mathbf{p}}(t,\mathbf{x}),$$

(5)

Then, let us focus on Equation (3) in the sector in which heavy quarks with a momentum of the order $\sqrt{MT}$ interact with gluons with an energy of the order T. This would be equivalent to performing the matching between NRQCD and ${\mathrm{NRQCD}}_{sh}$ at tree level

$${\mathcal{L}}_{\xi }=\sum _{\mathbf{p}\ne 0}{\xi }_{\mathbf{p}}^{†}\left\{i{D}_0-\frac{{(\mathbf{p}-i\nabla +g\mathbf{A})}^2}{2M}+c_4\frac{p^4}{8M^3}+c_{F}g\frac{\sigma ·\mathbf{B}}{2M}\right\}{\xi }_{\mathbf{p}}+\cdots$$

(6)

It is more convenient to rearrange the terms in the following way:

$${\mathcal{L}}_{\xi }=\sum _{\mathbf{p}\ne 0}{\xi }_{\mathbf{p}}^{†}\left\{i{D}_0-\frac{p^2}{2M}+i\frac{\mathbf{p}·(\nabla +ig\mathbf{A})}{M}+\frac{{(\nabla +ig\mathbf{A})}^2}{2M}+c_4\frac{p^4}{8M^3}+c_{F}g\frac{\sigma ·\mathbf{B}}{2M}\right\}{\xi }_{\mathbf{p}}+\cdots$$

(7)

The advantage of this equation compared to the previous one is that now each term has a well-defined power counting. Since the only scale that has not been integrated out is T, it follows that ${\xi }_{\mathbf{p}}$ is of size $T^{3/2}$. Therefore, the first two terms in Equation (7) are of the order 1. The third term is of the order $\sqrt{\frac{T}{M}}$ and the rest of the terms are of the order $\frac{T}{M}$. We have obtained Equation (7) by tree-level manipulations of the NRQCD Lagrangian. However, it is easy to convince ourselves that this corresponds to the heavy quark sector of ${\mathrm{NRQCD}}_{sh}$. The only differences that a proper matching would bring up in this case are possible sub-leading corrections to M, and the Wilson coefficients $c_4$ and $c_F$. Note that Galilean symmetry and reparametrization invariance [37] constrain the form of Equation (7).

Let us now mention some features of ${\mathrm{NRQCD}}_{sh}$ we believe are worth emphasizing. The first one is that, up to a trivial shift in the energy, the field ${\xi }_{\mathbf{p}}$ behaves at the leading-order as a static quark. This implies that we can use ${\mathrm{NRQCD}}_{sh}$ to compute properties of a heavy quark with a momentum of order $\sqrt{MT}$ by performing perturbations around the static case. The second remarkable feature that we want to comment on is that the study of the Wigner distribution is very much simplified in ${\mathrm{NRQCD}}_{sh}$. The Wigner distribution of a heavy quark in NRQCD takes the following form:

$$f_W(\mathbf{p},\mathbf{R})=\int d^{3}r{\mathrm{e}}^{-i\mathbf{p}\mathbf{r}}\mathit{Tr}\left({\psi }^{†}\left(t,\mathbf{R}-\frac{\mathbf{r}}{2}\right)\varphi \left(t:\mathbf{R}-\frac{\mathbf{r}}{2},\mathbf{R}+\frac{\mathbf{r}}{2}\right)\psi \left(t,\mathbf{R}+\frac{\mathbf{r}}{2}\right)\rho \right),$$

(8)

where $\varphi$ is a Wilson line introduced in the definition such that the Wigner distribution is gauge invariant [38] and $\rho$ is the density matrix. By virtue of this Wilson line, the $\mathbf{p}$ momentum in the Wigner distribution corresponds to the kinetic momentum and not to the canonical momentum conjugate [38]. We are interested in the case in which $p\sim \sqrt{MT}$ and the dependence of $f_W$ with R is only sizeable over distances much larger than $1/T$. In this case, $f_W$ is matched in ${\mathrm{NRQCD}}_{sh}$ into

$$\begin{array}{c}\hfill f_W(\mathbf{r},\mathbf{R})=Z_0\left(\mathbf{r}\right)\mathit{Tr}\left({\xi }_{\mathbf{p}}^{†}(t,\mathbf{R}){\xi }_{\mathbf{p}}(t,\mathbf{R})\rho \right)\\ \hfill +i{Z}_1\left(\mathbf{r}\right)\mathit{Tr}\left(\mathbf{r}·\mathbf{D}{\xi }_{\mathbf{p}}^{†}(t,\mathbf{R}){\xi }_{\mathbf{p}}(t,\mathbf{R})-\mathbf{r}·{\xi }_{\mathbf{p}}^{†}(t,\mathbf{R})\mathbf{D}{\xi }_{\mathbf{p}}(t,\mathbf{R})\rho \right)\\ \hfill +Z_2\left(\mathbf{r}\right)r\mathit{Tr}\left({\xi }_{\mathbf{p}}^{†}(t,\mathbf{R})\mathbf{r}·\mathbf{E}(t,\mathbf{R}){\xi }_{\mathbf{p}}(t,\mathbf{R})\rho \right)+\cdots \end{array}$$

(9)

where $f_W(\mathbf{r},\mathbf{R})$ is the Fourier transform of Equation (8). This equation is valid up to terms smaller than $\frac{T}{M}$. Note that, since both the lhs and the rhs can be understood as a pseudo-probability distribution normalized to 1, $Z_0\left(\mathbf{r}\right)=1+\mathcal{O}\left(r^2{T}^2\right)$. In conclusion, up to sub-leading corrections, we can identify $\mathit{Tr}\left({\xi }_{\mathbf{p}}^{†}(t,\mathbf{R}){\xi }_{\mathbf{p}}(t,\mathbf{R})\rho \right)$ with the Wigner distribution. The use of ${\mathrm{NRQCD}}_{sh}$ allows for seeing the Wigner distribution as a probability distribution (at least at the leading order) encoded in a local operator.

Let us now discuss some properties of the propagators in ${\mathrm{NRQCD}}_{sh}$ in thermal field theory that will be useful in the following. From now on, and in order to simplify the notation, we will drop the sub-index $\mathbf{p}$ from the field $\xi$. First, let us discuss the dilute limit, which we define as the limit in which

$$-\mathit{Tr}\left(\xi (t,{\mathbf{r}}_1){\xi }^{†}(t,{\mathbf{r}}_2)\rho \right)\gg \mathit{Tr}\left({\xi }^{†}(t,{\mathbf{r}}_2)\xi (t,{\mathbf{r}}_1)\rho \right).$$

(10)

In this limit,

$$S_{11}(t,\mathbf{r})=\mathit{Tr}\left(\mathcal{P}{\xi }_1(t,\mathbf{r}){\xi }_1^{†}(0,\mathbf{0})\rho \right)\sim \theta \left(t\right)\mathit{Tr}\left(\xi (t,\mathbf{r}){\xi }^{†}(0,\mathbf{0})\rho \right),$$

(11)

$$S_{22}(t,\mathbf{r})=\mathit{Tr}\left(\mathcal{P}{\xi }_2(t,\mathbf{r}){\xi }_2^{†}(0,\mathbf{0})\rho \right)\sim \theta (-t)\mathit{Tr}\left(\xi (t,\mathbf{r}){\xi }^{†}(0,\mathbf{0})\rho \right),$$

(12)

$$S_{12}(t,\mathbf{r})=S^{<}(t,\mathbf{r})=\mathit{Tr}\left(\mathcal{P}{\xi }_1(t,\mathbf{r}){\xi }_2^{†}(0,\mathbf{0})\rho \right)=-\mathit{Tr}\left({\xi }^{†}(0,\mathbf{0})\xi (t,\mathbf{r})\right)\rho )\sim 0,$$

(13)

and

$$S_{21}(t,\mathbf{r})=S^{>}(t,\mathbf{r})=\mathit{Tr}\left(\mathcal{P}{\xi }_2(t,\mathbf{r}){\xi }_1^{†}(0,\mathbf{0})\rho \right)\sim \mathit{Tr}\left(\xi (t,\mathbf{r}){\xi }^{†}(0,\mathbf{0})\rho \right).$$

(14)

We note that these results are true as long as we are exactly in the dilute limit. It is also useful to write the expression in the Keldysh representation [39].

$$S_R(t,\mathbf{r})=\theta \left(t\right)\mathit{Tr}\left(\left\{\xi (t,\mathbf{r}){\xi }^{†}(0,\mathbf{0})\right\}\rho \right)\sim S_{11}(t,\mathbf{r}),$$

(15)

$$S_A(t,\mathbf{r})=-\theta (-t)\mathit{Tr}\left(\left\{\xi (t,\mathbf{r}){\xi }^{†}(0,\mathbf{0})\right\}\rho \right)\sim -S_{22}(t,\mathbf{r}),$$

(16)

and

$$S_S(t,\mathbf{r})=\mathit{Tr}\left(\left[\xi (t,\mathbf{r}){\xi }^{†}(0,\mathbf{0})\right]\rho \right)\sim S_{21}(t,\mathbf{r}).$$

(17)

The following useful relation is exactly fulfilled in the dilute limit

$$S_S(t,\mathbf{r})\sim S_R(t,\mathbf{r})-S_A(t,\mathbf{r}).$$

(18)

This implies that, in the dilute limit, the retarded propagator contains all the relevant information, since the advanced propagator is the complex conjugate of the retarded one.

Let us now consider that the density of heavy quarks is small but non-zero. In this case, $S^{<}$ is small but not zero. Since this propagator is directly related to the distribution of heavy quarks, finite density corrections in other propagators have a sub-leading impact on the mentioned distribution.

Finally, let us write the heavy quark propagator in ${\mathrm{NRQCD}}_{sh}$ at tree level. The tree-level propagators are more compactly written in the Keldysh representation

$$S_R^0(t,\mathbf{r})=\int \frac{d^{4}k}{{\left(2\pi \right)}^4}{\mathrm{e}}^{-i{k}_{0}t+i\mathbf{k}\mathbf{r}}\frac{i}{k_0-\frac{p^2}{2M}+iϵ}=\theta \left(t\right){\mathrm{e}}^{-\frac{i{p}^{2}t}{2M}}{\delta }^{\left(3\right)}\left(\mathbf{r}\right),$$

(19)

$$S_A^0(t,\mathbf{r})=\int \frac{d^{4}k}{{\left(2\pi \right)}^4}{\mathrm{e}}^{-i{k}_{0}t+i\mathbf{k}\mathbf{r}}\frac{i}{k_0-\frac{p^2}{2M}-iϵ}=-\theta (-t){\mathrm{e}}^{-\frac{i{p}^{2}t}{2M}}{\delta }^{\left(3\right)}\left(\mathbf{r}\right),$$

(20)

and

$$\begin{array}{c}\hfill S_S^0(t,\mathbf{r},\mathbf{R})=\mathit{Tr}\left(\left[{\xi }^0\left(t,\mathbf{R}+\frac{\mathbf{r}}{2}\right){{\xi }^0}^{†}\left(0,\mathbf{R}-\frac{\mathbf{r}}{2}\right)\right]\rho \right)\\ \hfill =\int \frac{d^{4}k}{{\left(2\pi \right)}^4}{\mathrm{e}}^{-i{k}_{0}t+i\mathbf{k}\mathbf{r}}2\pi \delta \left(k_0-\frac{p^2}{2M}\right)(1-2f^0(\mathbf{p}+\mathbf{k},\mathbf{R})),\end{array}$$

(21)

where the 0 super-index denotes tree-level quantities and we have assumed that the tree-level distribution function of heavy quarks $f^0$ is a very smooth function over distances of size $\mathbf{r}$. Note, also, that to be consistent with the ${\mathrm{NRQCD}}_{sh}$ power counting, the function $f^0$ (or the resummed equivalent that we will introduce later) must be expanded

$$f^0(\mathbf{p}+\mathbf{k},\mathbf{R})\sim f^0(\mathbf{p},\mathbf{R})+\mathbf{k}{\nabla }_{\mathbf{p}}{f}^0(\mathbf{p},\mathbf{R})+\frac{k^i{k}^j}{2}{\mathsf{\Delta }}_{\mathbf{p}}^{ij}{f}^0(\mathbf{p},\mathbf{R})+\cdots$$

(22)

This implies that

$$\begin{array}{c}\hfill S_S^0(t,\mathbf{r},\mathbf{R})\sim {\mathrm{e}}^{-\frac{i{p}^{2}t}{2M}}\left[(1-2f^0(\mathbf{p},\mathbf{R})){\delta }^{\left(3\right)}\left(\mathbf{r}\right)+2\nabla {\delta }^{\left(3\right)}\left(\mathbf{r}\right){\mathbf{\nabla }}_{\mathbf{p}}{f}^0(\mathbf{p},\mathbf{R})\right.\\ \hfill \left.-{\mathsf{\Delta }}^{ij}{\delta }^{\left(3\right)}\left(\mathbf{r}\right){\mathsf{\Delta }}_{\mathbf{p}}^{ij}{f}^0(\mathbf{p},\mathbf{R})+\cdots \right].\end{array}$$

(23)

Finally, let us mention that the Feynmann rules for ${\mathrm{NRQCD}}_{sh}$ can be found in Appendix A.

3. Langevin Effective Theory

In this section, we introduce Langevin Effective Theory (LET). It is obtained from ${\mathrm{NRQCD}}_{sh}$ after integrating out degrees of freedom with an energy of the order of the scale T. In this case, we are dealing with heavy quarks with a momentum of $\mathbf{p}+\mathbf{k}$ where the residual momentum k is now much smaller than the temperature. Note that we can always redefine p such that this relation is fulfilled. In order to construct LET, we follow the observations of [19]. The influence functional of LET can be written as ${\mathrm{e}}^{i{S}_{LET}}$, with

$$S_{LET}=S_1-S_2+\mathcal{I}+S_{HTL},$$

(24)

where $S_1$ ($S_2$) is the piece of the action that involves only heavy quark fields of type 1(2). $\mathcal{I}$ is a new type of contribution in which heavy quark fields of both type 1 and 2 appear. Finally, $S_{HTL}$ is the well-known Hard Thermal Loop action [40]. If the scale T did not induce any dissipative effects, then $\mathcal{I}$ would be zero and $S_2$ would be equal to $S_1$, just changing the fields of type 2 by fields of type 1.

The construction of LET is substantially simplified by applying the dilute limit. We know that, in this limit, we can ignore the doubling of the degrees of freedom when computing Green functions in which only heavy fields of type 1 (or 2) appear [7,8]. Moreover, in the dilute limit, and as long as the heavy particle is non-relativistic in the frame in which the medium is at rest, the symmetries of ${\mathrm{NRQCD}}_{sh}$ are not broken by the presence of the medium. Therefore, the form of $S_1$ is equal to the action of ${\mathrm{NRQCD}}_{sh}$. However, the Wilson coefficients are different. In fact, some of these Wilson coefficients can be complex. Then,

$$\begin{array}{c}\hfill S_1=\int d^{4}x{\xi }_1^{†}\left[(1-\delta Z)\left(i{D}_0-\frac{p^2}{2M}\right)-\delta E+i\frac{\mathsf{\Gamma }}{2}+i\frac{\mathbf{p}·(\nabla +ig\mathbf{A})}{M}+\frac{{(\nabla +ig\mathbf{A})}^2}{2M}\right.\\ \hfill \left.+{\tilde{c}}_4\frac{p^4}{8M^3}+{\tilde{c}}_{F}g\frac{\sigma ·\mathbf{B}}{2M}\right]{\xi }_1+\cdots \end{array}$$

(25)

where we allow $\delta Z$, $\delta E$, and $\mathsf{\Gamma }$ to be polynomials involving M, T and $p^2$. Note that $\delta Z$ can have both a real and an imaginary part. The real part of $\delta Z$ can be reabsorbed by a unitary transformation. Whether the imaginary part of $\delta Z$ can be reabsorbed by a transformation of the fields is beyond the scope of this work. We note that, since our main focus is the study of $\frac{df}{dt}$, we can ignore a non-zero value of $\delta Z$. The reason is that the wave-function renormalization is not a secular effect, meaning that its effects on the evolution on the distribution of heavy quarks does not become larger as we study longer times. This is in contrast with what happens to corrections to the mass and the decay width. For example, even if the decay width is small, it becomes a leading order effect at long enough times.

$S_2$ can be obtained from $S_1$ by changing the fields of type 1 to fields of type 2 and by making the complex conjugate of the Wilson coefficients. In our case, this means making the changes $\delta Z\to \delta Z^{*}$ and $i\frac{\mathsf{\Gamma }}{2}\to -i\frac{\mathsf{\Gamma }}{2}$. We note also that, since $\delta E$ and $\mathsf{\Gamma }$ are obtained by performing a matching computation to ${\mathrm{NRQCD}}_{sh}$, they are polynomials with the following structure:

$$\delta E={\alpha }_{0}T+{\alpha }_1\frac{T^2}{M}+{\alpha }_2\frac{p^{2}T}{M^2}+\mathcal{O}\left(\frac{T^3}{M^2}\right),$$

(26)

$$\mathsf{\Gamma }={\beta }_{0}T+{\beta }_1\frac{T^2}{M}+{\beta }_2\frac{p^{2}T}{M^2}+\mathcal{O}\left(\frac{T^3}{M^2}\right).$$

(27)

Regarding $\mathcal{I}$, it can be fixed by imposing the following conditions:

• $S_{LET}$ is equal to zero if fields of type 1 are equal to fields of type 2 [19].
• In the dilute limit, the propagator $S_{12}$ is zero.

Using this, we obtain

$$\mathcal{I}=\int d^{4}x{\xi }_2^{†}\left[2i\mathrm{Im}Z\left(i{D}_0-\frac{p^2}{2M}\right)-i\mathsf{\Gamma }-i\mathrm{Im}{\tilde{c}}_4\frac{p^4}{4M^3}-i\mathrm{Im}{\tilde{c}}_F\frac{\sigma ·\mathbf{B}}{M}\right]{\xi }_1+\cdots$$

(28)

Let us now discuss the issue of gauge invariance. The equations we wrote are only invariant regarding transformations in which the fields are modified in the same in both branches of the Schwinger-Keldysh contour. More specifically, there is an explicit invariance under the following type of infinitesimal gauge transformations:

$$\begin{array}{c}\hfill {\xi }_1\to {\xi }_1+ig\mathsf{\Lambda }{\xi }_1,\\ \hfill {\xi }_2\to {\xi }_2+ig\mathsf{\Lambda }{\xi }_2,\\ \hfill A_1^{\mu }\to A_1^{\mu }+{\partial }^{\mu }\mathsf{\Lambda }+ig[\mathsf{\Lambda },A_1^{\mu }],\\ \hfill A_2^{\mu }\to A_2^{\mu }+{\partial }^{\mu }\mathsf{\Lambda }+ig[\mathsf{\Lambda },A_2^{\mu }].\end{array}$$

(29)

Apparently, something has been missed in going from ${\mathrm{NRQCD}}_{sh}$ to LET. ${\mathrm{NRQCD}}_{sh}$ is invariant under transformations in which each branch of the Schwinger-Keldysh contour is independently transformed with a different ${\mathsf{\Lambda }}_i$. How to recover this more general invariance in EFTs, where terms mixing the two branches appear, was discussed in [20,21,22]; however, the solution found is not suitable for our case. We have proposed an alternative solution in Appendix B, where we discuss this issue and its solution in more detail. However, at the end of the day, this issue has little practical importance for the computation at hand and will only complicate the notation. Careful readers might have noticed that we also did not mention whether the gauge fields entering $D_0$ and $\mathbf{B}$ in Equation (28) are of type 1 or 2. More details about this are also given in Appendix B.

The previous way of presenting the action of $S_{LET}$ is useful in achieving matching to the full theory. However, we might rearrange the contributions in a more physically meaningful way.

$$S_{LET}={\tilde{S}}_1-{\tilde{S}}_2+\mathcal{D},$$

(30)

where

$$\begin{array}{c}\hfill {\tilde{S}}_i=\int d^{4}x{\xi }_i^{†}\left[(1-\mathrm{Re}\delta Z)\left(i{D}_0-\frac{p^2}{2M}\right)-\delta E+i\frac{\mathbf{p}·(\nabla +ig\mathbf{A})}{M}+\frac{{(\nabla +ig\mathbf{A})}^2}{2M}\right.\\ \hfill \left.+\mathrm{Re}{\tilde{c}}_4\frac{p^4}{8M^3}+\mathrm{Re}{\tilde{c}}_{F}g\frac{\sigma ·\mathbf{B}}{2M}\right]{\xi }_i+\cdots \end{array}$$

(31)

is the unitary part of the evolution and $\mathcal{D}$ is the dissipative part

$$\begin{array}{c}\hfill \mathcal{D}=i\int d^{4}x\left(\begin{array}{cc}{\xi }_1^{†}& {\xi }_2^{†}\end{array}\right)\left[-\mathrm{Im}\delta Z\left(i{D}_0-\frac{p^2}{2M}\right)+\frac{\mathsf{\Gamma }}{2}+\mathrm{Im}{\tilde{c}}_4\frac{p^4}{8M^3}+\mathrm{Im}{\tilde{c}}_F\frac{\sigma ·\mathbf{B}}{2M}\right]\\ \hfill \otimes \left(\begin{array}{cc}1& 0\\ -2& 1\end{array}\right)\left(\begin{array}{c}{\xi }_1\\ {\xi }_2\end{array}\right).\end{array}$$

(32)

It might be illustrative to write $S_{LET}$ in the Keldysh basis [39,41]. Let us introduce

$${\xi }_{Cl}=\frac{{\xi }_1+{\xi }_2}{2},$$

(33)

which we call the classical field and

$${\xi }_Q={\xi }_1-{\xi }_2,$$

(34)

the quantum field. There are two remarkable properties that can be seen just introducing this basis. First, the condition that $S_{LET}$ is zero when fields of type 1 are equal of type 2 can be rephrased as imposing that $S_{LET}=0$ when ${\xi }_Q=0$. It also follows that a QFT describing a closed system (without dissipation) only has terms containing an odd number of quantum fields. However, in an EFT obtained after integrating out medium degrees of freedom (such as the one we are studying in this paper), we might have terms containing an even number of quantum fields. In summary, terms with an even number of quantum fields are forbidden in ${\tilde{S}}_i$ while they are allowed in $\mathcal{D}$. To see this more explicitly, let us introduce the following operators:

$$\begin{array}{c}\hfill {\mathsf{\Delta }}_S=(1-\mathrm{Re}\delta Z)\left(i{D}_0-\frac{p^2}{2M}\right)-\delta E+i\frac{\mathbf{p}·(\nabla +ig\mathbf{A})}{M}+\frac{{(\nabla +ig\mathbf{A})}^2}{2M}\\ \hfill +\mathrm{Re}{\tilde{c}}_4\frac{p^4}{8M^3}+\mathrm{Re}{\tilde{c}}_{F}g\frac{\sigma ·\mathbf{B}}{2M},\end{array}$$

(35)

and

$${\mathsf{\Delta }}_D=-\mathrm{Im}\delta Z\left(i{D}_0-\frac{p^2}{2M}\right)+\frac{\mathsf{\Gamma }}{2}+\mathrm{Im}{\tilde{c}}_4\frac{p^4}{8M^3}+\mathrm{Im}{\tilde{c}}_F\frac{\sigma ·\mathbf{B}}{2M}.$$

(36)

Then, we can write $S_{LET}$ in the following way:

$$S_{LET}=\int d^{4}x\left({\xi }_{Cl}^{†}{\mathsf{\Delta }}_S{\xi }_Q+{\xi }_Q^{†}{\mathsf{\Delta }}_S{\xi }_{Cl}+{\xi }_Q^{†}{\mathsf{\Delta }}_D{\xi }_{Cl}-{\xi }_{Cl}^{†}{\mathsf{\Delta }}_D{\xi }_Q+{\xi }_Q^{†}{\mathsf{\Delta }}_D{\xi }_Q\right).$$

(37)

Writing $S_{LET}$ in this way provides some extra insight. In a situation in which ${\xi }_Q$ is suppressed, the leading contribution comes from terms linear in ${\xi }_Q$. Integrating over ${\xi }_Q$ considering only these linear terms gives a Dirac delta that forces ${\xi }_{Cl}$ to follow the classical equations of motion. Terms quadratic in ${\xi }_Q$ can be seen as originating from a classical random source; therefore, we can understand them as fluctuations [42].

Now, let us discuss the case in which we still consider that heavy quarks are dilute but we take into account the first corrections proportional to their density. We call this NLO dilute corrections. Regarding the symmetries of the EFT, we consider that $f\left(\mathbf{p}\right)$ does not have any preferred direction other than $\mathbf{p}$ itself. Taking this into account, $\mathcal{D}$ has to be modified in order to include an extra term

$$\begin{array}{c}\hfill \mathcal{D}=i\int d^{4}x\left(\begin{array}{cc}{\xi }_1^{†}& {\xi }_2^{†}\end{array}\right)\left[\left(-\mathrm{Im}\delta Z\left(i{D}_0-\frac{p^2}{2M}\right)+\frac{\mathsf{\Gamma }}{2}+\mathrm{Im}{\tilde{c}}_4\frac{p^4}{8M^3}+\mathrm{Im}{\tilde{c}}_F\frac{\sigma ·\mathbf{B}}{2M}\right)\right.\\ \hfill \left.\otimes \left(\begin{array}{cc}1& 0\\ -2& 1\end{array}\right)+\mathsf{\Delta }\mathsf{\Gamma }\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)\right]\left(\begin{array}{c}{\xi }_1\\ {\xi }_2\end{array}\right).\end{array}$$

(38)

We note that the term that we have added is the leading one that we can add that does not fulfill the condition that the $S_{12}$ propagator has to be zero (dilute limit) but that fulfills the rest of the conditions that we have discussed in this section. Again, we can rewrite the action in the classical-quantum basis,

$$\begin{array}{c}\hfill S_{LET}=\int d^{4}x\left({\xi }_{Cl}^{†}{\mathsf{\Delta }}_S{\xi }_Q+{\xi }_Q^{†}{\mathsf{\Delta }}_S{\xi }_{Cl}+{\xi }_Q^{†}{\mathsf{\Delta }}_D{\xi }_{Cl}-{\xi }_{Cl}^{†}{\mathsf{\Delta }}_D{\xi }_Q+{\xi }_Q^{†}{\mathsf{\Delta }}_D{\xi }_Q\right.\\ \hfill \left.+{\xi }_{Cl}^{†}\mathsf{\Delta }\mathsf{\Gamma }{\xi }_Q-{\xi }_Q^{†}\mathsf{\Delta }\mathsf{\Gamma }{\xi }_{Cl}\right).\end{array}$$

(39)

We note that all Wilson coefficients can be affected by the NLO dilute corrections in a sub-leading way. The specific property of $\mathsf{\Delta }\mathsf{\Gamma }$ is that it vanishes in the exact dilute limit.

There are additional symmetries in the EFT that impose relations between the different terms in $\mathsf{\Gamma }$ and $\mathsf{\Delta }\mathsf{\Gamma }$. The origin of these relations is the Schwinger-Keldysh symmetry [20,21,22], also known as the fluctuation–dissipation theorem. However, let us postpone this discussion until the next section, since it is very much related to the evolution of $f\left(\mathbf{p}\right)$.

4. Evolution of $\mathit{f}\left(\mathbf{p}\right)$ in LET

In this section, we show how to compute the evolution of $f\left(\mathbf{p}\right)$ in the EFT we introduced before. At the same time, this will also allow us to introduce further constraints on the Wilson coefficients by imposing that the fluctuation–dissipation theorem is fulfilled1. In more physical terms, this means the following. We are studying the case of a dilute distribution of heavy quarks evolving in a large bath in thermal equilibrium at a temperature such that $T\ll M$. Therefore, at very long times, $f\left(\mathbf{p}\right)$ must be equal to the thermal distribution, ${\mathrm{e}}^{-\frac{p^2}{2MT}}$.

The information on the distribution of heavy quarks can be found more directly in the $S^{<}\left(p\right)$ propagator, which, as we have seen before, goes to zero in the dilute limit. In order to be more precise, we define the distribution function f, such that

$$S^{<}\left(k\right)=f(\mathbf{p}+\mathbf{k})(S_R\left(k\right)-S_A\left(k\right)),$$

(40)

where, as usual, we are refering to the propagators of the field ${\xi }_p$. To study this propagator, it is convenient to use the Kadanoff–Baym equations [43]. A recent application of these equations in the context of heavy-ion collisions can be found in [44]. These equations are deduced by performing a Dyson–Schwinger type of resummation of the self-energies. In our case, we will just perform a resummation of the tree-level self-energies as obtained directly from the LET influence functional.

$$\begin{array}{c}\hfill {\mathsf{\Pi }}_R=\delta E-i\frac{\mathsf{\Gamma }}{2},\\ \hfill {\mathsf{\Pi }}_A=\delta E+i\frac{\mathsf{\Gamma }}{2},\\ \hfill {\mathsf{\Pi }}^{<}=-i\mathsf{\Delta }\mathsf{\Gamma }.\end{array}$$

(41)

The first and second lines of the previous equation follow directly from Equation (37). The third line follows from Equation (38). We note that, when studying a system evolving in a plasma, some kind of resummation is always needed to deal with secular effects. In other words, small perturbations that grow with time need to be resummed because they become leading order effects at long enough times.

The propagator $S^{<}$ is a function of two times,

$$S^{<}(t_1,{\mathbf{r}}_1;t_2,{\mathbf{r}}_2)=\mathit{Tr}\left({\xi }_1(t_1,{\mathbf{r}}_1){\xi }^{†}(t_2,{\mathbf{r}}_2)\rho \right).$$

(42)

In thermal equilibrium and due to translational invariance, it is only a function of $\tau =t_1-t_2$. More generally, it is also a function of $t=\frac{t_1+t_2}{2}$. In order to study the evolution of $f\left(\mathbf{p}\right)$, the more direct way is to look at the evolution of $S^{<}$ as a function of t for $\tau =0$. The Kadanoff–Baym equations give us the evolution on each time separately

$$\begin{array}{c}\hfill {\partial }_{t_1}{S}^{<}=-i\left(\frac{p^2}{2M}+{\mathsf{\Pi }}_R\right)S^{<}+i{\mathsf{\Pi }}^{<}{S}_A,\\ \hfill {\partial }_{t_2}{S}^{<}=i\left(\frac{p^2}{2M}+{\mathsf{\Pi }}_A\right)S^{<}-i{\mathsf{\Pi }}^{<}{S}_R.\end{array}$$

(43)

From this, it follows that the evolution with t is given by

$${\partial }_t{S}^{<}=-i({\mathsf{\Pi }}_R-{\mathsf{\Pi }}_A)S^{<}-i{\mathsf{\Pi }}^{<}(S_R-S_A).$$

(44)

Note that the spectral function $\rho$ is given by $\rho =S_R-S_A$, which, at tree level, is ${\rho }_0=2\pi \delta \left(k_0-\frac{p^2}{2M}\right)$. We are studying the case of a heavy particle interacting with a medium in thermal equilibrium. At very long times, we should arrive to a steady state in which the distribution of heavy particles is also in thermal equilibrium. This means that at late times

$$S^{<}=-\frac{{\mathsf{\Pi }}^{<}(S_R-S_A)}{{\mathsf{\Pi }}_R-{\mathsf{\Pi }}_A}.$$

(45)

At the same time, at thermal equilibrium, the fluctuation–dissipation theorem must be fulfilled

$$S^{<}=f_{eq}(\mathbf{p}+\mathbf{k})(S_R-S_A),$$

(46)

where $f_{eq}\left(p\right)=N{\mathrm{e}}^{-\frac{p^2}{2M}}$ is the thermal equilibrium distribution function and N is a normalization factor. We obtain the relation

$${\mathsf{\Pi }}^{<}(\mathbf{p},\mathbf{k},f_{eq})=f_{eq}(\mathbf{p}+\mathbf{k})({\mathsf{\Pi }}_R(\mathbf{p},\mathbf{k})-{\mathsf{\Pi }}_A(\mathbf{p},\mathbf{k})),$$

(47)

that must be fulfilled. Note that we have explicitly written the dependency with the momentum and the distribution function. To simplify, let us focus on the case $\mathbf{k}=\mathbf{0}$ and write the expression in terms of the Wilson coefficients of LET

$$\mathsf{\Delta }\mathsf{\Gamma }\left(\mathbf{p},{\mathrm{e}}^{-\frac{p^2}{2M}}\right)={\mathrm{e}}^{-\frac{p^2}{2M}}\mathsf{\Gamma }\left(\mathbf{p}\right).$$

(48)

$\mathsf{\Delta }\mathsf{\Gamma }$ is proportional to f or its derivatives

$$\mathsf{\Delta }\mathsf{\Gamma }(\mathbf{p},f_{eq})=\left({\gamma }^0+\frac{{\gamma }^1}{M}+\frac{{\gamma }^2{p}^2}{M^2}\right)f_{eq}\left(\mathbf{p}\right)+\frac{{\gamma }^3}{M}\mathbf{p}{\nabla }_{\mathbf{p}}\left(f_{eq}\left(\mathbf{p}\right)\right)+{\gamma }^4{\mathsf{\Delta }}_{\mathbf{p}}{f}_{eq}\left(\mathbf{p}\right)+\mathcal{O}{\left(\frac{T}{M}\right)}^{3/2}.$$

(49)

At the same time, a similar expansion is valid for $\mathsf{\Gamma }$

$$\mathsf{\Gamma }={\beta }_0+\frac{{\beta }_1}{M}+\frac{{\beta }_2{p}^2}{M^2}.$$

(50)

Then, this implies that the following relations must be fulfilled:

$$\begin{array}{c}\hfill {\gamma }_0={\beta }_0,\\ \hfill {\gamma }_1-\frac{3{\gamma }_4}{T}={\beta }_1,\\ \hfill {\gamma }_2-\frac{{\gamma }_3}{T}+\frac{{\gamma }_4}{T^2}={\beta }_2.\end{array}$$

(51)

In order to compute $\frac{\partial f\left(\mathbf{p}\right)}{\partial t}$, we can use Equations (44), (46) and (51) in the case $\mathbf{k}=\mathbf{0}$

$$\frac{\partial f\left(\mathbf{p}\right)}{\partial t}=-\mathsf{\Gamma }f\left(\mathbf{p}\right)+\mathsf{\Delta }\mathsf{\Gamma }f\left(\mathbf{p}\right)=\left(\frac{3}{M}-\frac{p^2}{M^{2}T}\right)\left(\frac{{\gamma }_4}{T}-{\gamma }_3\right)f\left(\mathbf{p}\right)+\frac{{\gamma }_3}{M}\nabla \left(\mathbf{p}f\left(\mathbf{p}\right)\right)+{\gamma }_4{\mathsf{\Delta }}_{\mathbf{p}}f\left(\mathbf{p}\right).$$

(52)

Another condition that must be fulfilled is unitarity. This implies that

$$\frac{\partial }{\partial t}\int d^{3}pf\left(\mathbf{p}\right)=0.$$

(53)

This imposes an additional condition, this is ${\gamma }_4=T{\gamma }_3$. Note that what we obtain is the evolution of $f\left(\mathbf{p}\right)$ that corresponds to a Fokker–Planck equation, in which the drag coefficient is given by $\frac{{\gamma }_3}{M}$.

In summary, imposing all the constraints coming from unitarity and the fluctuation–dissipation theorem, we obtain the following expression for $\mathsf{\Delta }\mathsf{\Gamma }$

$$\mathsf{\Delta }\mathsf{\Gamma }(\mathbf{p},f)=\left({\beta }_0+\frac{{\beta }_1}{M}+\frac{{\beta }_2{p}^2}{M^2}\right)f\left(\mathbf{p}\right)+\frac{\kappa }{2MT}{\nabla }_{\mathbf{p}}\left(\mathbf{p}f\left(\mathbf{p}\right)\right)+\frac{\kappa }{2}{\mathsf{\Delta }}_{\mathbf{p}}f\left(\mathbf{p}\right)+\mathcal{O}{\left(\frac{T}{M}\right)}^{3/2},$$

(54)

where $\kappa$ is the heavy quark diffusion coefficient [45]. We see that the matching between ${\mathrm{NRQCD}}_{sh}$ and $LET$ becomes substantially simplified. To determine $\mathsf{\Delta }\mathsf{\Gamma }$ at the order we are interested in, we only need to know $\mathsf{\Gamma }$ in the dilute limit and the value of the heavy quark diffusion coefficient. Moreover, if we are only interested in the evolution of $f\left(\mathbf{p}\right)$, we only need to know $\kappa$, as expected, and we can see by writing the evolution of $f\left(\mathbf{p}\right)$, considering all the constraints, that

$$\frac{\partial f\left(\mathbf{p}\right)}{\partial t}=\frac{\kappa }{2MT}\nabla \left(\mathbf{p}f\left(\mathbf{p}\right)\right)+\frac{\kappa }{2}{\mathsf{\Delta }}_{\mathbf{p}}f\left(\mathbf{p}\right).$$

(55)

5. Matching between ${\mathbf{NRQCD}}_{\mathit{sh}}$ and LET

In this section, we are going to perform matching between ${\mathrm{NRQCD}}_{sh}$ and LET in the one-gluon exchange approximation. The reader might wonder why we do not compute the matching at a given order in perturbation theory instead. The reason is related to the special properties of gauge theories at a finite temperature. A one-loop matching will lead to $\kappa =0$, meaning that $f\left(\mathbf{p}\right)$ does not change. On the other hand, a two-loop matching would involve a complex computation of $\mathsf{\Gamma }$ that has little phenomenological impact due to the symmetries that we discussed in the previous section. The one-gluon exchange approximation has the virtue of giving a finite contribution to all Wilson coefficients already at the leading order, and we can regard it as an intermediate step between a one-loop and a two-loop matching.

The strategy to perform the matching is going to be the following. First, we are going to compute ${\mathsf{\Pi }}_R$ in ${\mathrm{NRQCD}}_{sh}$; from this we can obtain $\delta M$ and $\mathsf{\Gamma }$. Second, we are going to compute ${\mathsf{\Pi }}^{<}$, but only the piece proportional to ${\nabla }_{\mathbf{p}}f\left(\mathbf{p}\right)$, since we have seen in Section 4 that this is enough to fix all the relevant Wilson coefficients. Regarding the computation of ${\mathsf{\Pi }}_R$, we are going to use the Keldysh basis following the graphical notation of [46]2. For the computation of ${\mathsf{\Pi }}^{<}$, we find it more convenient to use the $1-2$ basis. We note that, in this section, we are going to use the Feynmann rules discussed in Appendix A.

Let us begin by computing the leading contribution to ${\mathsf{\Pi }}_R$ in ${\mathrm{NRQCD}}_{sh}$; in other words, the terms that in the EFT power counting are of the order T. In Figure 2, we show the relevant diagrams. If we write $\delta E$ as

$$\delta E={\alpha }_0+\frac{{\alpha }_1}{M}+\frac{{\alpha }_2{p}^2}{M^2},$$

(56)

and we use the notation of Equation (50), then the contribution from the diagrams in Figure 2 is

$$\begin{array}{c}\hfill {\alpha }_0=-\frac{i{g}^2{C}_F}{4}\int \frac{d^{3}q}{{\left(2\pi \right)}^3}({\mathsf{\Delta }}_{00,R}+{\mathsf{\Delta }}_{00,A})(0,\mathbf{q}),\\ \hfill {\beta }_0=2T{g}^2{C}_F\int \frac{d^{3}q}{{\left(2\pi \right)}^3}{\rho }_{00}^{\prime }(0,\mathbf{q}),\end{array}$$

(57)

where ${\mathsf{\Delta }}_{00}$ is the temporal gluon propagator and ${\rho }^{\prime }$ is the first derivative over the temporal component of the spectral function of the temporal gluon propagator [16]. Our computation is performed in the Coulomb gauge.

Next, let us discuss the possible corrections to $\delta E$ and $\mathsf{\Gamma }$ of the order $T\sqrt{\frac{T}{M}}$. By symmetry arguments, there should not be any. However, let us check explicitly that this is the case. The relevant diagrams are shown in Figure 3. Each circled correction to the propagator introduces a factor $\frac{\mathbf{p}\mathbf{q}}{M}$, where $\mathbf{q}$ is the spatial momentum of the temporal gluon. We can take $\mathbf{p}$ out of the integrand, and since the medium does not have any preferred direction in space, the result of the integration can only be zero.

Regarding the diagrams contributing at the order $\frac{T^2}{M}$, we can divide them into two classes: those that contribute to ${\alpha }_1$ and ${\beta }_1$, and those contributing to ${\alpha }_2$ and ${\beta }_2$. The first class is shown in Figure 4. They give the following contributions:

$$\begin{array}{c}\hfill {\alpha }_1=\frac{g^2{C}_F}{4}\int \frac{d^{4}q}{{\left(2\pi \right)}^4}{\mathsf{\Delta }}_{ii,S}\left(q\right)+\frac{g^2{C}_F}{8}\int \frac{d^{4}q}{{\left(2\pi \right)}^4}{q}^2\left(\frac{1}{{(q_0-iϵ)}^2}+\frac{1}{{(q_0+iϵ)}^2}\right){\mathsf{\Delta }}_{00,S}\left(q\right)\\ \hfill -\frac{i{g}^2{C}_F}{8}\int \frac{d^{3}q}{{\left(2\pi \right)}^3}{q}^2\frac{d({\mathsf{\Delta }}_{00}^R+{\mathsf{\Delta }}_{00}^A)}{d{q}_0}(0,\mathbf{q}),\\ \hfill {\beta }_1=-\frac{g^2{C}_f}{4}\int \frac{d^{3}q}{{\left(2\pi \right)}^3}{q}^2{\rho }_{00}^{\prime }(0,\mathbf{q}).\end{array}$$

(58)

Diagrams contributing to ${\alpha }_2$ and ${\beta }_2$ are shown in Figure 5. They result in the following contributions:

$$\begin{array}{c}\hfill {\alpha }^2=-i\frac{g^2{C}_F}{12}\int \frac{d^{3}q}{{\left(2\pi \right)}^3}({\mathsf{\Delta }}_{ii,R}+{\mathsf{\Delta }}_{ii,A}),\\ \hfill {\beta }^2=\frac{g^2{C}_{F}T}{6}\int \frac{d^{3}q}{{\left(2\pi \right)}^3}{\rho }_{ii}^{\prime }(0,\mathbf{q})-\frac{g^2}{6}\int \frac{d^{4}q}{{\left(2\pi \right)}^4}{\left(\frac{i}{q_0+iϵ}\right)}^3{\mathsf{\Delta }}_{00,S}\left(q\right){\mathbf{q}}^2.\end{array}$$

(59)

Let us now compute ${\mathsf{\Pi }}^{<}$ in ${\mathrm{NRQCD}}_{sh}$ in order to perform the matching. We are going to focus on the piece proportional to $\mathbf{p}{\nabla }_{\mathbf{p}}f\left(\mathbf{p}\right)$ since, as we discussed previously, it is the only extra term needed to perform the matching. To obtain this, we need to go to the second order in the expansion shown in Equation (23). Moreover, we have to take from the ${\mathrm{NRQCD}}_{sh}$ Lagrangian the interaction similar to $\frac{{\xi }^{†}\mathbf{p}\nabla \xi }{M}$. Therefore, we need to focus on the diagram shown in Figure 6. Note that, for the computation of ${\mathsf{\Pi }}^{<}$, we are using the $1-2$ basis instead of the Schwinger-Keldysh one. From this matching, we obtain

$$\kappa =\frac{g^2{C}_{F}T}{3}\int \frac{d^{3}q}{{\left(2\pi \right)}^3}{q}^2{\rho }_{00}^{\prime }(0,\mathbf{q}).$$

(60)

We can check that this result is compatible with the perturbative computation of $\kappa$ in the Coulomb gauge [25,45,47]. Note that, in the $q_0\to 0$ limit, we can approximate $\mathbf{E}$ by $-\nabla A_0$, as long as we are not using the unitary gauge. It is also important to take into account that we are enganging in a slight abuse of language by identifying the transport coefficient $\kappa$ with the Wilson coefficient of LET. Strictly speaking, the Wilson coefficient only includes the contribution to $\kappa$ from the scale we are integrating out, in this case T.

With this, we have finished the matching to the order of interest. However, for the sake of cross-checking some of the relations we discussed in the previous section, we are going to discuss an alternative way to match the value of $\kappa$. This can be achieved looking at the diagram in Figure 7. However, now we have to take the third term in Equation (23) to get the contribution proportional to $\mathsf{\Delta }f\left(\mathbf{p}\right)$. Doing this, we can check that we obtain exactly the same value of $\kappa$. There is even a third way of computing $\kappa$ which is to look at the terms in $\mathsf{\Delta }\mathsf{\Gamma }$ that function similarly to $\frac{f\left(\mathbf{p}\right)}{M}$ and subtract ${\beta }^1$. However, we will not discuss further this way of achieving the matching of $\kappa$.

Now we are ready to discuss how to derive a Fokker–Planck equation for heavy quarks in QCD including the effects of the scales T and $gT$. The first step has already been discussed. We match ${\mathrm{NRQCD}}_{sh}$ to LET to include the effect of the scale T in the Wilson coefficients. The second step is to compute ${\mathsf{\Pi }}_R$ and ${\mathsf{\Pi }}^{<}$ in LET, including the effects of the scale $gT$. Note, however, that the same symmetries and cancellations apply now, and, therefore, we only need to compute the contribution of the scale $gT$ to $\kappa$. This can be obtained by computing the diagram in Figure 6, but now the loop has to be computed using the HTL gluon propagator. Proceeding in this way, we obtain the perturbative estimate of $\kappa$. This corresponds to Equation (60) but now using the two loop perturbative result for $q\sim T$ and the HTL approximation for $q\sim m_D$ [25,47].

6. Heavy Quarks with Momentum $\mathit{P}\sim \mathit{T}$

In this section, we discuss the case of a heavy quark with a tri-momentum of order T. In this case, we will proceed as follows:

• As a starting point, we will use NRQCD Lagrangian.
• We will integrate out the scale T. After this, we will arrive again to LET. However, now the Wilson coefficients and the power counting will be slightly different, since p is now of the order T instead of $\sqrt{MT}$ as before.
• In this case, $\mathsf{\Delta }\mathsf{\Gamma }$ will have a more involved dependency with p. As a result, instead of a Fokker–Planck equation, we obtain an evolution similar to a Boltzmann equation.

Let us start looking at LET. The Lagrangian is again formed by combining Equations (25) and (28) in the same way as before. The main difference is the value of the Wilson coefficients.

Regarding $\delta E$ and $\mathsf{\Gamma }$, it is easy to check that the results are exactly the same as in the case $p\sim \sqrt{MT}$. The reason is that the diagrams that appear and carry out the respective matchings in ${\mathrm{NRQCD}}_{sh}$ (for $p\sim \sqrt{MT}$) and in NRQCD (for $p\sim T$) are the same at any order in the coupling constant, the only difference being the relative size of each term. However, the situation is different for $\mathsf{\Delta }\mathsf{\Gamma }$; the reason for this is that now this Wilson coefficient is, in general, a non-trivial function of f and p and we cannot assume a polynomial form like in Equation (54). However, we can assume that it is a polynomial in inverse powers of $1/M$ from the structure of NRQCD. Then, we can write

$$\mathsf{\Delta }\mathsf{\Gamma }(\mathbf{p},f)=\int \frac{d^{3}q}{{\left(2\pi \right)}^3}\left(C_1\left(q\right)+\frac{C_2\left(q\right)q^2}{M}+\frac{C_3\left(q\right)(2\mathbf{q}\mathbf{p}-q^2)}{M}\right)f(\mathbf{p}-\mathbf{q}).$$

(61)

There are several constraints that can be imposed on $C_i$. Equation (48) must be fulfilled, but now we need to take into account that we need to expand ${\mathrm{e}}^{-\frac{p^2}{2MT}}$ to the order desired since $\frac{p^2}{2MT}\ll 1$. This implies that

$$\begin{array}{c}\hfill \int \frac{d^{3}q}{{\left(2\pi \right)}^3}{C}_1\left(q\right)={\beta }_0,\\ \hfill \int \frac{d^{3}q}{{\left(2\pi \right)}^3}{q}^2\left(C_2\left(q\right)-C_3\left(q\right)-\frac{C_1\left(q\right)}{2T}\right)={\beta }_1,\end{array}$$

(62)

Additional constraints come from unitarity. We can repeat the arguments of Section 4 to arrive at

$${\partial }_{t}f\left(\mathbf{p}\right)=\mathsf{\Delta }\mathsf{\Gamma }(\mathbf{p},f)-\mathsf{\Gamma }f\left(\mathbf{p}\right).$$

(63)

Imposing that ${\partial }_t\int \frac{d^{3}p}{{\left(2\pi \right)}^3}f(\mathbf{p}=0)$, we deduce that

$$\int \frac{d^{3}p}{{\left(2\pi \right)}^3}\mathsf{\Delta }\mathsf{\Gamma }(\mathbf{p},f)=\mathsf{\Gamma }.$$

(64)

Starting from the previous equation and using the change of variables $\mathbf{p}\to \mathbf{p}+\mathbf{q}$, we obtain the additional constraint

$$\int \frac{d^{3}q}{{\left(2\pi \right)}^3}{q}^2\left(C_2\left(q\right)+C_3\left(q\right)\right)={\beta }_1.$$

(65)

Using the previous equation and some trivial manipulations, we arrive at our final form of the evolution equation

$$\begin{array}{c}\hfill {\partial }_{t}f\left(\mathbf{p}\right)=\int \frac{d^{3}q}{{\left(2\pi \right)}^3}\left(C_1\left(q\right)+\frac{C_2\left(q\right)q^2}{M}\right)(f(\mathbf{p}-\mathbf{q})-f\left(\mathbf{p}\right))\\ \hfill +\frac{1}{M}\int \frac{d^{3}q}{{\left(2\pi \right)}^3}{C}_3\left(q\right)(2\mathbf{q}\mathbf{p}-q^2)(f(\mathbf{p}-\mathbf{q})+f\left(\mathbf{p}\right)).\end{array}$$

(66)

As an illustration, let us discuss the perturbative matching. This can be carried out by looking at the diagrams in Figure 6 and Figure 7, but now taking into account that the expansion in Equation (23) cannot be achieved since $p\sim q$. We obtain that

$$\begin{array}{c}\hfill C_1\left(q\right)=g^2{C}_{F}T{\rho }_{00}^{\prime }(0,\mathbf{q}),\\ \hfill C_3\left(q\right)=-\frac{g^2{C}_F{\rho }_{00}^{\prime }(0,\mathbf{q})}{4},\\ \hfill C_2\left(q\right)=0.\end{array}$$

(67)

In fact, the last equation seems to be true at any perturbative order just because of the structure of NRQCD. We can see explicitly that the conditions in Equations (62) and (65) are fulfilled.

Finally, let us discuss the physical picture that the EFT treatment of the cases $p\sim T$ and $p\sim \sqrt{MT}$ provides when they are put together. If we start with a heavy quark at rest with the medium, it will start to gain momentum due to collisions. At this stage, the evolution of the system is given by a Boltzmann equation like Equation (66). This evolution is such that there is almost an equal probability to gain momentum than to lose it. Therefore, we will naturally arrive at a situation in which $p\gg T$ after some time. In this case, the evolution can already be described by a Langevin equation. As the momentum keeps increasing, the drag force starts to become important. Finally, we arrive to an equilibrium distribution in which momenta much larger than $\sqrt{MT}$ are very rare. Let us note that the arguments needed to arrive at these qualitative statements are valid, even if the medium is strongly coupled.

7. Discussion and Conclusions

In this manuscript, we have discussed an EFT approach to the derivation of the evolution of a heavy quark in a medium. Our main motivation was to pave the way for future developments in the study of quarkonium suppression. We would like to obtain a Lindblad equation for quarkonium without assuming a weakly coupled plasma and using the hierarchy of energy scales that appears in the problem. This has already been achieved in the case that the medium sees quarkonium as a small color dipole [10,11]. However, we would like to generalize to the case $Tr\sim 1$, as this would allow modeling more realistically excited states of bottomonium and charmonium at LHC energies. This is challenging because the EFT in which heavy quarkonium is best described is pNRQCD, and it is obtained after integrating out the scale $1/r$. However, if $Tr\sim 1$, the matching between NRQCD and pNRQCD is modified by the medium. The appearance of an imaginary part of the potential indicates that the medium induces dissipative effects that cannot be encoded in an effective action; instead, we need an influence functional in which terms mixing the two branches of the Schwinger-Keldysh contour appear. Therefore, we need to modify pNRQCD to include this kind of term. The motivation of this manuscript was to study a simpler case in order to pave the way for this modification of pNRQCD in the future.

The simpler case we studied is a single heavy quark interacting with a medium. We achieved this following the EFT philosophy, in which each energy scale must be treated separately using a series of EFTs. We distinguished two cases: the case $p\sim \sqrt{MT}$ and the case $p\sim T$. In both cases, we start with NRQCD; however, in the first case, we integrate out the scale $\sqrt{MT}$ as an intermediate step going from NRQCD to ${\mathrm{NRQCD}}_{sh}$. Finally, we integrate the scale T in both cases to go to a new EFT that we named LET. The structure of this EFT is given by Equation (39), which is one of the main results of this manuscript. This formula combines the constraints of the Keldysh symmetry together with the special features of heavy particles in the dilute limit. Further constraints are provided by the fluctuation–dissipation theorem, which, together with the EFT power counting, leads to Equation (55) for $p\sim \sqrt{MT}$ and to Equation (66) for $p\sim T$. This is the second main result from this manuscript: that Equations (55) and (66) can be obtained simply by symmetry and power counting arguments within the EFT approach. The physical picture is also quite transparent; a heavy quark with a momentum of the order T can be described by a Boltzmann equation and will slowly increase its momentum until we reach a moment in which the momentum is large enough so that a Langevin equation becomes valid.

Let us now discuss in more detail the relevance of the study presented in this manuscript in the context of quarkonium suppression. We have seen that the state of a heavy quark cannot be described by a single position x; instead, we need the position in the upper (lower) Schwinger-Keldysh branch $x_1$, ($x_2$). Note that $x_1-x_2$ is related to p by the Wigner transform. The increase of p we discussed below is seen in coordinate space as a narrowing of the $x_1-x_2$ coordinate or, put differently, the density matrix becoming almost diagonal in the coordinate basis. This behavior has also been observed in quarkonium [48,49]. One main difference between the case of quarkonium and that of a single heavy quark is that the number of coordinates needed to describe the system doubles. In quarkonium, we can talk about the center-of-mass coordinate, R, and the relative coordinate, r. However, in non-equilibrium situations, we need to consider the doubling of degrees of freedom, and then we have to take into account four coordinates, $R_1$, $R_2$, $r_1$, and $r_2$. A complete EFT treatment of quarkonium would identify all the possible scale hierarchies that can be constructed with these coordinates and that appear in the problem of quarkonium suppression. The case studied in this manuscript corresponds to a situation in which the quark and the antiquark are very far apart, and we have studied the cases $x_1-x_2\sim \frac{x_1+x_2}{2}$ and $x_1-x_2\ll \frac{x_1+x_2}{2}$. In [10,11], the case studied corresponds to $r_1-r_2\sim \frac{r_1+r_2}{2}\ll \frac{1}{T}$ and $T\gg E$. Perturbative studies have shown that the Lindblad equation that appears in the case $T\gg E$ naturally leads to a decrease of $r_1-r_2$. This produces a change in the power counting that, among other things, makes the drag force not a perturbation (NLO corrections in $E/T$ that include the drag force were added in [50,51]). Regarding the relation between the scales $1/r$, T and E, pNRQCD has been studied for all possible relations [7,8,9] but only taking into account the part analogous to $S_i$ in the notation of Equation (24). This is enough to discuss spectroscopy, but not to study the evolution of quarkonium inside a plasma. What is missing is the construction of the analogous to $\mathcal{I}$ in pNRQCD, something that we hope this work has helped to develop.

Finally, let us discuss possible extensions of the present work in the context of heavy quark propagation itself. As in any EFT, an obvious improvement would be to compute higher order corrections in the $\sqrt{\frac{T}{M}}$ expansion to improve our knowledge about the evolution of $f\left(\mathbf{p}\right)$. Another possible direction is to relax the assumptions about the symmetries of the problem. For example, until now, we have assumed that both the medium and the distribution of heavy quarks are homogeneous in space and isotropic. It would be interesting to relax these conditions, as they are not completely fulfilled in heavy-ion collisions. Finally, the case $p\gg \sqrt{MT}$ is also interesting. In this case, we would be studying the case of a heavy quark that loses energy until it thermalizes with the medium.