An Effective Field Theory study of medium heavy quark evolution

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.


I. INTRODUCTION
The 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, the 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 is the contribution 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 ∼ Mv with v ≪ 1, and the binding energy E ∼ Mv 2 .Non-relativistic QCD (NRQCD) [2,3] is an EFT valid for energies smaller than M. Since M ≫ Λ 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 done, we can use these Wilson coefficients to compute the properties of quarkonium states with finite heavy quark masses.This is a rigorous way to show 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 has 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 ≫ T [10][11][12][13].Regarding jets, there are recently developed EFTs for the study of jet broadening in [14] and jet sub-structure 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 finite temperature, computations of non-static properties are done using the so-called realtime formalism [16].In this formalism, there is a doubling of 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 fig. 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 φ.Then the action can be written as where the subindex C means that the integration is done 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 have access to the distribution function we would need to compute φ(t)φ(0) , where the ordering of the operators is as written.Note also that, although in eq. ( 1) there are no terms containing at the same time both type of fields, the propagator of the field φ i is not diagonal in the 1, 2 indices [16].
What we have explained until now regarding the Schwinger-Keldysh contour (for the specific case of a QFT describing the evolution of a field φ) is standard textbook material.
However, it was important to remind it 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 eq. ( 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 can not be encoded in a Hamiltonian [17].In other words, in general it does not exist an operator 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 the path integral formalism.
The effect of high energy fields on low energy ones can not be encoded in an effective action when we are dealing with an open quantum system.We need a more general object called influence functional instead [18].For the example of a φ field where P means ordering along the Schwinger-Keldysh contour and F is the influence functional.Only in the case of an isolated system it happens that , 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 and EFT.Up to 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-relativisty systems.In this case, the effects of the doubling of degrees of freedom are neutralized in the approximation in which the heavy particles are dilute [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 order of √ T M , being T the temperature.
• The temperature T and other energy scales that the medium might induced.
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 these steps: 1.As a starting point, we can use NRQCD.This automatically encodes the effects of degrees of freedom with energy of order M.
2. It is straight-forward to define an EFT for energies smaller than √ 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 NRQCD for semi-hard fields (NRQCD sh ) [34].This EFT was originally introduced to study modes with energy of the order MΛ QCD in the study of heavy quarkonium at T = 0.In our case, the role of Λ QCD is played by T .However, symmetries and power counting arguments remain the same.
3. 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.
4. We perform the matching between NRQCD sh and LET in the one-gluon exchange approximation as an illustration.This allows to encode the influence of degrees of freedom with energy of order T .
5. 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 ∼ 1 r .The problem that we study in this manuscript is simpler than quarkonium in the regime T ∼ 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 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 finite temperature can be also done 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 ∼ 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 ≫ T .However, the Wilson coefficients are different.We observe that there is a smooth transition between the cases p ∼ T and p ≫ T , as it should.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 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 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 II, we will discuss NRQCD sh .Next, in section III , we will present LET.Section IV discusses the computation of the evolution of the distribution of heavy quarks within LET.In section V, we discuss the matching between NRQCD sh and LET.In the section VI, we consider the case of a heavy quark with momentum p ∼ T .Finally, in section VII, we give our conclusions.

II. NRQCD sh
Our starting point is the NRQCD Lagrangian [2,3].To fix the notation, we write its heavy quark sector Here D µ = ∂ µ − igA µ , σ 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 √ MT plus a possible residual momentum of order T .The gauge field is only sensitive to the scale T .In eq. ( 3) we have only written terms of order ψ † ψ T 2 M or lower.These are the terms that we need in order to compute 1 f df dt to order T 2 M where f (p, 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 order √ MT while each interaction with the medium changes the momentum of the heavy quark by an amount of order T .This implies that the term ψ † D 2 2M ψ 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 T M or 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 it is done in SCET and OSEFT.Let us divide the spatial momentum of heavy quarks in two pieces where p is of order √ MT and k is a residual momentum of order T .To improve this we can perform the following transformation Then, let us focus on eq.(3) in the sector in which heavy quarks with momentum of order √ MT interact with gluons with energy of order T .This would be equivalent to performing the matching between NRQCD and NRQCD sh at tree level It is more convenient to rearrange the terms in the following way 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 ξ p is of size T 3/2 .Therefore, the first two terms in eq. ( 7) are of order 1.The third term is of order T M and the rest of terms are of order T M .We have obtained eq. ( 7) by tree-level manipulations of the NRQCD Lagrangian.However, it is easy to convince our selves that it corresponds to the heavy quark sector of 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] constraint the form of eq. ( 7).
Let us now mention some features of NRQCD sh we believe are worth emphasizing.The first one is that, up to a trivial shift in the energy, the field ξ p behaves at leading-order as a static quark.This implies that we can use NRQCD sh to compute properties of a heavy quark with a momentum of order √ MT by performing perturbations around the static case.
The second remarkable feature that we want to comment is that the study of the Wigner distribution is very much simplified in NRQCD sh .The Wigner distribution of a heavy quark in NRQCD takes the following form where φ is a Wilson line introduced in the definition such that the Wigner distribution is gauge invariant [38] and ρ is the density matrix.By virtue of this Wilson line, the 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 ∼ √ MT and the dependence of f W with R is only sizeable over distances much larger than 1/T .In this where f W (r, R) is the Fourier transform of eq. ( 8).This equation is valid up to terms smaller than T M .Note that, since both the lhs and the rhs can be understood as a pseudoprobability distribution normalized to 1, Z 0 (r) = 1 + O(r 2 T 2 ).In conclusion, up to subleading corrections, we can identify Tr ξ † p (t, R)ξ p (t, R)ρ with the Wigner distribution.The use of NRQCD sh allows seeing the Wigner distribution as a probability distribution (at least at leading order) encoded in a local operator.
Let us now discuss some properties of the propagators in NRQCD sh in thermal field theory that will be useful in the following.Now on, and in order to simplify the notation, we will drop the sub-index p from the field ξ.First, let us discuss the dilute limit, which we define as the limit in which In this limit, and 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]. and The following useful relation is exactly fulfilled in the dilute limit 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.
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 with 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 NRQCD sh at tree level.The tree level propagators are more compactly written in the Keldysh representation and 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 r.Note also that to be consistent with the NRQCD sh power counting, the function f 0 (or the resummed equivalent that we will introduce later) must be expanded This implies that Finally, let us mention that the Feynmann rules for NRQCD sh can be found in Appendix A.

III. LANGEVIN EFFECTIVE THEORY
In this section we introduce Langevin Effective Theory (LET).It is obtained from 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 momentum p + k where now the residual momentum k is 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 e iS LET , with where S 1 (S 2 ) is the piece of the action that involves only heavy quark fields of type 1 (2).
I is a new type of contribution in which heavy quark fields of both type 1 and 2 appear.
Finally, S HT L is the well-known Hard Thermal Loop action [40].If the scale T did not induce any dissipative effects, then 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 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 NRQCD sh are not broken by the presence of the medium.Therefore, the form of S 1 is equal to the action of NRQCD sh .However, the Wilson coefficients are different.In fact, some of these Wilson coefficients can be complex.Then, where we allow both δZ, δE and Γ to be polynomials involving M, T and p 2 .Note that δZ can have both a real and an imaginary part.The real part of δZ can be reabsorbed by a unitary transformation.Whether the imaginary part of δ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 df dt , we can ignore a non-zero value of δZ.The reason is that the wavefunction 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 to 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 large 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 δZ → δZ * and i Γ 2 → −i Γ 2 .We note also, that since δE and Γ are obtained by performing a matching computation to NRQCD sh , they are polynomials with the following structure Regarding 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 get 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: Apparently, something has been missed in going from NRQCD sh to LET.NRQCD sh is invariant under transformations in which each branch of the Schwinger-Keldysh contour is independently transformed with a different Λ 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 the 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 notived that we also did not mention whether the gauge fields entering D 0 and B in eq. ( 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 to perform the matching to the full theory.However, we might rearrange the contributions in a more physically meaningful way. where is the unitary part of the evolution and D is the dissipative part It might be illustrative to write S LET in the Keldysh basis [39,41].Let us introduce which we call the classical field and 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 ξ 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 Si while they are allowed in D. To see this more explicitly, let us introduce the following operators: and Then, we can write S LET in the following way Writing S LET in this way provides some extra insight.In a situation in which ξ Q is suppressed, the leading contribution comes from terms linear in ξ Q .Integrating over ξ Q considering only these linear terms gives a Dirac delta that forces ξ Cl to follow the classical equations of motion.Terms quadratic in ξ 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 (p) does not have any preferred direction other than p itself.Taking this into account, D has to be modified in order to include an extra term 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 conditions that we have discussed in this section.Again, we can rewrite the action in the classical -quantum basis, We note that all Wilson coefficients can be affected by the NLO dilute corrections in a sub-leading way.The specific property of ∆Γ is that it vanishes in the exact dilute limit.
There are additional symmetries in the EFT that impose relations between the different terms in Γ and ∆Γ.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 next section, since it is very much related with the evolution of f (p).

IV. EVOLUTION OF f (p) IN LET
In this section, we show how to compute the evolution of f (p) 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 fulfilled [43].On 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 ≪ M. Therefore, at very large times f (p) must be equal to the thermal distribution, e − p 2 2M T .
The information on the distribution of heavy quarks can be found more directly in the S < (p) 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 where, as usual, we are refering to the propagators of the field ξ p .To study this propagator, it is convenient to use the Kadanoff-Baym equations [44].A recent application of these equations in the context of heavy-ion collisions can be found in [45].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.
The first and second lines of the previous equation follow directly from eq. ( 37).The third line follows from eq. ( 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 large enough times.
The propagator S < is a function of two times, In thermal equilibrium and due to translational invariance it is only a function of τ = t 1 −t 2 .
More generally, it is also a function of t = t 1 +t 2 2 .In order to study the evolution of f (p), the more direct way is to look at the evolution of S < as a function of t for τ = 0.The Kadanoff-Baym equations give us the evolution on each time separately From this, it follows that the evolution with t is given by Note that the spectral function ρ is given by ρ = S R − S A , which at tree level is ρ 0 = 2πδ k 0 − p 2 2M .We are studying the case of a heavy particle interacting with a medium in thermal equilibrium.At very large 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 At the same time, at thermal equilibrium the fluctuation-dissipation theorem must be fulfilled where f eq (p) = Ne − p 2 2M is the thermal equilibrium distribution function and N a normalization factor.We obtain the relation that must be fulfilled.Note that we have explicitly written the dependency with the mo- ∆Γ is proportional to f or its derivatives At the same time, a similar expansion is valid for Γ Then, this implies that the following relations must be fulfilled: In order to compute ∂f (p) ∂t , we can use eqs.( 44), ( 46) and ( 51) Another condition that must be fulfilled is unitarity.This implies that This imposes an additional condition, this is γ 4 = T γ 3 .Note that what we obtain is the evolution of f (p) that corresponds to a Fokker-Planck equation, in which the drag coefficient is given by γ 3 M .In summary, imposing all the constraints coming from unitarity and the fluctuationdissipation theorem, we get the following expression for ∆Γ where κ is the heavy quark diffusion coefficient [46].We see that the matching between NRQCD sh and LET gets substantially simplified.To determine ∆Γ at the order we are interested, we only need to know Γ in the dilute limit and the value of the heavy quark diffusion coefficient.Moreover, if we are only interested in the evolution of f (p) we only need to know κ, as it was expected, and we can see by writing the evolution of f (p) considering all the constraints ∂f (p) V. MATCHING BETWEEN NRQCD sh AND LET In this section, we are going to perform the matching between 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 finite temperature.A one-loop matching will lead to κ = 0, meaning that f (p) does not change.On the other hand, a two-loop matching would involve a complex computation of Γ 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 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 Π R in NRQCD sh , from this we can obtain δM and Γ.Second, we are going to compute Π < , but only the piece proportional to ∇ p f (p), since we have seen in section IV that this is enough to fix all the relevant Wilson coefficients.Regarding the computation of Π R , we are going to use the Keldysh basis following the graphical notation of [47][48].For the computation of Π < , 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 Π R in NRQCD sh , in other words, the terms that in the EFT power counting are of order T .In fig. 2 we show the relevant diagrams.If we write δE as and we use the notation of eq. ( 50), then the contribution from the diagrams in fig. 2 is where ∆ 00 is the temporal gluon propagator and ρ ′ 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 possible corrections to δE and Γ of order T T M .By symmetry arguments, there should not be any.However, let us check explicitly that it is the case.The relevant diagrams are shown in fig. 3.Each circled correction to the propagator introduces a factor pq M , where q is the spatial momentum of the temporal gluon.We can take 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 order T 2 M , we can divide them in two classes.Those that contribute to α 1 and β 1 , and those contributing to α 2 and β 2 .The first class is shown in fig. 4. They give the following contribution: Diagrams contributing to α 2 and β 2 are shown in fig. 5.They result in the following contribution: Let us now compute Π < in NRQCD sh in order to perform the matching.We are going to focus on the piece proportional to p∇ p f (p) since, as we discussed previously, it is the only extra term needed to perform the matching.To obtain this, we need to go to second order in the expansion shown in eq.(23).Moreover, we have to take from the NRQCD sh Lagrangian the interaction that goes like ξ † p∇ξ M .Therefore, we need to focus on the diagram shown in fig.6.Note that for the computation of Π < we are using the 1 − 2 basis instead of the Schwinger-Keldysh one.From this matching, we obtain We can check that this result is compatible with the perturbative computation of κ in the Coulomb gauge [25,46,49].Note that in the q 0 → 0 limit we can approximate E by −∇A 0 , as long as we are not using the unitary gauge.It is also important to take into account that we are making a small abuse of language by identifying the transport coefficient κ with the Wilson coefficient of LET.Strictly speaking, the Wilson coefficient only includes the contribution to κ 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 κ.This can be done looking at the diagram in fig. 7.However, now we have to take the third term in eq. ( 23) to get the contribution proportional to ∆f (p).Doing this, we can check that we get exactly the same value of κ.
There is even a third way of computing κ which is to look at the terms in ∆Γ that go like M and subtract β 1 .However, we are not going to discuss further this way of doing the matching of κ.
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 NRQCD sh to LET to include the effect of the scale T in the Wilson coefficients.
The second step is to compute Π R and Π < 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 κ.This can be obtained by computing the diagram in fig.6, but now the loop has to be computed using the HTL gluon propagator.
Proceeding in this way, we obtain the perturbative estimate of κ.It corresponds to eq. ( 60) but now using the two loop perturbative result for q ∼ T and the HTL approximation for q ∼ m D [25,49].

VI. HEAVY QUARKS WITH MOMENTUM p ∼ 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 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 now p is of order T instead of √ MT as before.
• In this case, ∆Γ will have a more involved dependency with p.As a result, instead of a Fokker-Planck equation, we get an evolution similar to a Boltzmann equation.
Let us start looking at LET.The Lagrangian is again formed by combining eq. ( 25) and eq.( 28) in the same way as before.The main difference is the value of the Wilson coefficients.
Regarding δE and Γ, it is easy to check that the results are exactly the same as in the case p ∼ √ MT .The reason is that the diagrams that appear doing the respective matchings in NRQCD sh (for p ∼ √ MT ) and in NRQCD (for p ∼ 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 ∆Γ, the reason is that now this Wilson coefficient is, in general, a non-trivial function of f and p and we can not assume a polynomial form like in eq. ( 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 There are several constraints that can be imposed to C i .Eq. ( 48) must be fulfilled, but now we need to take into account that we need to expand e − p 2 2M T to the order desired since 2M T ≪ 1.This implies that Additional constraints come from unitarity.We can repeat the arguments of section IV to arrive to Imposing that ∂ t d 3 p (2π) 3 f (p = 0) we deduce that Starting from the previous equation and using the change of variables p → p + q we obtain the additional constraint Using the previous equation and some trivial manipulations, we arrive at our final form of the evolution equation As an illustration, let us discuss the perturbative matching.This can be done by looking at the diagrams in figs.7 and 6, but now taking into account that the expansion in eq.(23) can not be done since p ∼ q.We obtain that Actually, 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 eqs.( 62) and (65) are fulfilled.
Finally, let us discuss the physical picture that the EFT treatment of the cases p ∼ T and p ∼ √ 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 eq. ( 66).This evolution is such that there is almost an equal probability to gain momentum than to lose it.Therefore, we will naturally arrive to a situation in which p ≫ 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 √ MT are very rare.Let us note that the arguments needed to arrive to these qualitative statements are valid even if the medium is strongly coupled.

VII. 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 done in the case that the medium sees quarkonium as a small color dipole [10,11].However, we would like to generalize to the case T r ∼ 1 as this would allow modeling more realistically excited states of bottomonium and charmonium at LHC energies.This is challenging because the EFT in 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 nonequilibrium 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 ∼ x 1 +x 2 2 and x 1 − x 2 ≪ x 1 +x 2 2 .In [10,11], the case studied corresponds to r 1 − r 2 ∼ r 1 +r 2 2 ≪ 1 T and T ≫ E. Perturbative studies have shown that the Lindblad equation that appears in the case T ≫ 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 [52,53]).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 analogous to S i part in the notation of eq. ( 24).This is enough to discuss spectroscopy, but not to study the evolution of quarkonium inside a plasma.What is missing, and we hope this work helps to develop, is the construction of the analogous to I in pNRQCD.
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 T M expansion to improve our knowledge about the evolution of f (p).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 ≫ √ MT is also interesting.In this case, we would be studying the case of a heavy quark that loses energy until thermalizing with the medium.

FIG. 1 .
FIG. 1. Representation of the Schwinger-Keldysh contour.The thicker line represents the complex time contour along which the path integral is defined in non-equilibrium field theory.It goes from −∞ to ∞ and then goes back to −∞ but decreasing the time's imaginary part by an infinitely small amount.

2 2M= e − p 2 2M
mentum and the distribution function.To simplify, let us focus on the case k = 0 and write the expression in terms of the Wilson coefficients of LET ∆Γ p, e − p Γ(p) .

FIG. 2 .
FIG.2.Diagrams contribution to Π R at order T in NRQCD sh .The computation is done in the Keldysh basis.

FIG. 3 .
FIG.3.Diagrams contribution to Π R at order T T M in NRQCD sh .

ACKNOWLEDGMENTS
MAE wants to thanks the careful reading of a first version of this manuscript and the useful comments of Juan Torres-Rincon, Joan Soto, Nora Brambilla and Antonio Vairo.The work of MAE has been supported by the Maria de Maetzu excellence program under project CEX2019-000918-M, by the Spanish Research State Agency under project PID2019-105614GB-C21 and by the grant 2021-SGR-249 of Generalitat de Catalunya.