Real Time Quarkonium Transport Coefficients in Open Quantum Systems from Euclidean QCD

Recent open quantum system studies showed that quarkonium time evolution inside the quark-gluon plasma is determined by transport coefficients that are defined in terms of a gauge invariant correlator of two chromoelectric field operators connected by an adjoint Wilson line. We study the Euclidean version of the correlator for quarkonium evolution and discuss the extraction of the transport coefficients from this Euclidean correlator, highlighting its difference from other problems that also require reconstructing a spectral function, such as the calculation of the heavy quark diffusion coefficient. Along the way, we explain why the transport coefficient $\gamma_{\rm adj}$ differs from $\gamma_{\rm fund}$ at finite temperature at $\mathcal{O}(g^4)$, in spite of the fact that their corresponding spectral functions differ only by a temperature-independent term at the same order. We then discuss how to evaluate the Euclidean correlator via lattice QCD methods, with a focus on reducing the uncertainty caused by infrared renormalons in determining the renormalization factor nonperturbatively.


I. INTRODUCTION
The scientific mission of relativistic heavy ion colliders is to investigate properties of the deconfined phase of nuclear matter in the high temperature regime, known as the quark-gluon plasma (QGP).In current heavy ion collision experiments, the QGP only lives for a short time period (roughly 10 fm/c in the laboratory frame) and we cannot directly measure its properties.Therefore, we use probes such as particle multiplicities and azimuthal distributions, jets and hadrons containing heavy quarks to indirectly study its properties.Various properties of the QGP are encoded in terms of gauge invariant correlation functions of field operators that often define transport coefficients showing up in the time evolution equations of the probes in the medium.Well-known examples include the shear viscosity (defined as a correlator of stressenergy tensors), the jet quenching parameter (a correlator of light-like Wilson lines) and the heavy quark diffusion coefficient (a correlator of two chromoelectric fields dressed with Wilson lines).Since the QGP is a strongly coupled fluid, nonperturbative determinations of these transport coefficients are crucial in our understanding of the QGP and QCD at finite temperature.Common nonperturbative methods include lattice QCD calculations and the holographic correspondence [1].One can also extract these transport coefficients from experimental data by solving in-medium evolution equations (which can be model dependent) for different values of the transport coefficients and then performing a Bayesian analysis [2][3][4][5][6].
Recently, thanks to the advance in applying the open quantum system framework to study jets [7] and quarkonia [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26] in the QGP (for recent reviews, see Refs. [27][28][29][30]), a novel correlator of two chromoelectric fields dressed with Wilson lines that determines transport properties of quarkonium in the medium was constructed [12,20].This correlator for quarkonium transport is similar to but different from the correlator defining the heavy quark diffusion coefficient [31,32] in terms of the ordering of the fields contained in the Wilson lines.Perturbative calculations in R ξ gauge showed that the spectral function of the correlator for quarkonium transport [33] differs from that for heavy quark transport [34] by a temperature independent constant at nextto-leading order (NLO).However, if both calculations had been performed in temporal axial gauge (A 0 = 0), one would, at first sight, have concluded that the two correlators were identical.This resulted in a puzzle: Since both correlators are defined in a gauge invariant way, calculations with different gauge choices must give the same result.This puzzle was resolved in Ref. [35], establishing the difference between the two correlators on a more solid ground in QCD.Beyond NLO, the heavy quark diffusion coefficient has been studied by using hard-thermalloop resummation [36], as well as nonperturbatively via the lattice QCD method [37][38][39][40][41] and the AdS/CFT correspondence [31,42,43].On the other hand, a recent AdS/CFT calculation showed that the analog quarkonium transport coefficients in N = 4 supersymmetric Yang-Mills (SYM) theory are zero [44], in stark contrast to the heavy quark diffusion coefficient value of √ λπT 3 at large coupling λ = g 2 N c ≫ 1.This difference is surprising because the heavy quark and quarkonium transport coefficients are defined by similar chromoelectric field correlators.Therefore, it is well motivated to study the quarkonium transport properties nonperturbatively in QCD.It is also crucial and urgent, since quarkonium production serves as an important probe of the QGP that is produced strongly coupled in current heavy ion collision experiments.
In this article, we discuss how to extract the quarkonium transport coefficients from lattice QCD calculations of a specific Euclidean chromoelectric correlator.The paper is organized as follows: We will first review the quarkonium transport coefficients in the real-time formalism in Section II, which are defined in terms of a correlator of two chromoelectric fields connected via an adjoint Wilson line.Then, in Section III we will discuss the Euclidean version of the correlator and how to relate it to its real time counterpart.Next, in Section IV the setup of a lattice QCD calculation of this Euclidean correlator will be discussed, with a focus on how to renormalize it.Finally, we will conclude and present our outlook in Section V.

II. QUARKONIUM TRANSPORT PROPERTIES
The quarkonium transport coefficients are defined in terms of time-ordered chromoelectric field operators, dressed with Wilson lines [12]: where ⟨O⟩ T ≡ Tr(O e −βH )/Tr(e −βH ), E a i is a chromoelectric field, W ab (t, 0) denotes a time-like Wilson line in the adjoint representation from t = 0 to t, T represents the time-ordering symbol, N c = 3 is the number of colors, and T F = 1/2 is the normalization of the fundamental representation generator matrices.To simplify the notation we have neglected the spatial coordinates, which are the same for all the fields, and will do so throughout the paper, unless the spatial coordinates are no longer the same.Both κ adj and γ adj appear in the Lindblad equation describing the time evolution of a heavy quarkantiquark pair (Q Q) at a small distance in the quantum Brownian motion limit [11,12]: where ρ S is the subsystem density matrix of the Q Q pair, γ adj ∆h S is the thermal correction to the vacuum Q Q Hamiltonian H S and L αi denotes the relevant Lindblad "jump" operators.Their explicit expressions are given in Appendix A. The κ adj parameter in the non-Hermitian part of the Lindblad equation determines the rate of transition between a Q Q pair in the color singlet state and that in the color octet state, as well as wavefunction decoherence.On the other hand, the γ adj parameter in the Hermitian part of the Lindblad equation controls the modification of the Q Q potential caused by the medium.
One way to interpret the integrations in Eq. ( 1) is as Fourier transforms that convert the time domain to the frequency domain.Consequently, the coefficients κ adj and γ adj are the zero frequency limits of frequencydependent correlation functions.Moreover, their behavior at finite frequency also turns out to be physically important.To explain the physical meaning of these correlation functions at finite frequency, we introduce pathordered chromoelectric field correlators [20,33] and consider their Fourier transforms The path-ordered version is more convenient to use at finite frequency and is consistent with the time-ordered version: It has been shown that κ adj = [g ++ adj ] > (ω = 0) [35].We note that because of the explicit operator ordering, only in [g ++ adj ] > the adjoint Wilson lines can be rewritten as W ab (t, 0), which appears in the time-ordered correlator shown in Eq. ( 1).Furthermore, one can obtain the time-ordered correlator that enters the definition of κ adj and γ adj by considering The path-ordered correlators at finite frequency appear in the Boltzmann (rate) equation for quarkonium dissociation and recombination, which is derived in the quantum optical limit of the open quantum system approach [20,33]: where n b (t, x) is the density of the quarkonium state b at time t, Γ denotes the dissociation rate and F represents the formation of the quarkonium state b from a recombining pair of unbound heavy quarks where a nonzero energy difference between the bound and unbound states ∆E = p2 rel /M + |E b | determines how the finite frequency dependence of the correlators appears in the transition rates.Here M is the heavy quark mass and E b is the binding energy of the quarkonium state b.The transition occurs via a color dipole interaction ⟨ψ b |r|Ψ p rel ⟩ between a bound Q Q state |ψ b ⟩ and an unbound scattering wave |Ψ p rel ⟩. f Q Q denotes the distribution of unbound heavy quark pairs with center-ofmass positions x and momenta p cm and relative positions x rel = 0 and momenta p rel .
The chromoelectric field correlator for quarkonium transport is different from that for heavy quark diffusion.In particular, the heavy quark diffusion coefficient κ fund and an analogous quantity γ fund (whose physical meaning has not been explored for heavy quark transport) are defined by where E i = E a i T a F is the Lie algebra-valued chromoelectric field, with the fundamental representation generator matrices normalized as Tr c (T a F T b F ) = T F δ ab .Also, Tr c denotes trace over color indices and U (t, 0) represents a time-like fundamental Wilson line from t = 0 to t.The subscript T in the expectation value denotes that the state on which this expectation value is calculated is a thermal density matrix, while the subscript Q means that this thermal state contains a static external color charge in the fundamental representation, e.g., a heavy quark.Mathematically, this expectation value is defined as and thus is different from that in Eq. (1).The fundamental Wilson line along the imaginary time at t = −∞ indicates the inclusion of the heavy quark effect on the thermal density matrix of the whole system.It is noted that the operators involved in the definition of κ fund and γ fund are path-ordered.We want to emphasize that the crucial difference between Eqs. (1) and ( 8) is not the representations of the Wilson lines, but the different orderings of the operators.

III. EUCLIDEAN CORRELATORS AND TRANSPORT COEFFICIENTS
As is well known, lattice QCD methods can only calculate correlation functions in Euclidean space and thus cannot be applied directly to study the real-time correlators defined in Eq. (3).In this section, we will introduce a Euclidean version of the correlator for quarkonium transport and discuss how to extract the quarkonium transport coefficients from the evaluation of this Euclidean correlator.As we will show, both the Euclidean correlator itself and the method to extract the quarkonium transport coefficients are different from the case of heavy quark diffusion in subtle and important aspects.To make the comparison more explicit, and also to take advantage of the apparent similarities between them, we will first review the extraction of the heavy quark diffusion coefficient from the corresponding Euclidean correlator.

A. Heavy Quark Diffusion
The Euclidean correlator relevant for the heavy quark diffusion case is given by [32] where β = 1/T is the inverse of the QGP temperature and ⟨•⟩ T = Tr(•e −βH )/Tr(e −βH ), with H the Hamiltonian of the QGP in the absence of any external color source.It has been shown that the heavy quark transport coefficient can be obtained from G fund (τ ) via [32,45] where the spectral function ρ fund (ω) is related to the Euclidean correlator through1 This correlator is constructed such that the standard Kubo-Martin-Schwinger (KMS) and analytic continuation relations hold as in textbook thermal field theory.Given an analytic expression for Gfund (ω n ), with ω n = 2πT n, n ∈ Z the Matsubara frequencies, one can extract the spectral function by taking the real part 2 of the retarded correlator obtained by analytic continuation ω n → −i(ω + iϵ) of this Euclidean correlator.This has been done both at weak [34] (QCD) and strong [31] (N = 4 SYM) coupling.However, at physical values of the coupling in QCD, the only tool available at the moment is lattice gauge theory, and as such, the reconstruction of the spectral function ρ fund through the relation (11) has received much attention in recent years [41,46,47].
Comparatively, the theoretical treatment of quarkonium transport coefficients has received less attention.We now aim to fill in this gap, and subsequently, to provide a recipe to determine these transport coefficients from lattice QCD calculations.To this end, we need to first construct a Euclidean version of the correlator for quarkonium transport that can be calculated via lattice QCD methods, and then explain how to extract the quarkonium transport coefficients from the evaluation of such an Euclidean correlator.We will answer these two questions in the following two subsections.Details of the lattice calculation of the Euclidean correlator will be discussed in the next section.

B. Euclidean Correlator for Quarkonium Transport
To construct the Euclidean correlator for quarkonium transport, we first note that because of the operator ordering in the definitions (3), we can equivalently write To perform the analytic continuation, it is best to explicitly isolate the t dependence from the field operators and write it purely in terms of time evolution factors.We let H be the Hamiltonian of the thermal bath QGP in the absence of any external color charge.When an external adjoint color charge is present, the Hamiltonian of the thermal bath is given by [H1 − gA c 0 (0)T c adj ] ab .The reason for the appearance of this modified Hamiltonian can be seen from converting the adjoint Wilson line back to the Schrödinger picture from the interaction picture Eq. ( 13) has the following physical interpretation: during the time interval between 0 and t the QGP evolves in the presence of an adjoint color charge, which is manifest in the modification of the Hamiltonian by −gA 0 .It is essentially a local modification to Gauss's law 3 , revealing the presence of a color octet Q Q pair.Outside this time interval the QGP evolves in the absence of external color sources.Using Eq. ( 13), one can write: = Tr H [e −βH ] , 3 An interesting question one can ask of this expression is whether we still have explicit gauge invariance.The answer is, naturally, affirmative.However, this is not as easy to see when considering time-dependent gauge transformations as it is for timeindependent gauge transformations.This is because the Hamiltonian also changes if one considers time-dependent gauge transformations, which is something to keep in mind when quantizing the theory.We will not pursue this further here, and we shall assume that H is already determined.For a thorough discussion on the quantization of gauge theories, we refer the reader to Ref. [48].
where the trace Tr H runs over physical states of the QGP.The analytic continuation is now direct, because all of the time dependence is in the exponentials.We just set t → −iτ , and identify the Euclidean gauge field A 4 with the Minkowski one by A 0 (0) = iA 4 (0) (which in turn means that the electric field picks up a factor of i), to find Tr H e Hτ E a i (0) e −(H−gA c where P denotes path-ordering.That is to say, we have proven that one of the real-time correlations we want to evaluate is related to an Euclidean correlation function by [g ++ adj ] > (−iτ ) = G adj (τ ).We note that the absence of the denominator term as in Eq. ( 9) is a result of the absence of a Wilson line along the imaginary time direction at t = −∞ in the definition of [g ++ adj ] > .In quarkonium dissociation, the initial state is a color singlet, whereas in heavy quark diffusion, the initial state is in a color triplet representation, whose effect appears explicitly in the initial thermal state.

C. Extraction of Quarkonium Transport
Coefficients from Euclidean QCD Now we discuss how to extract the quarkonium transport coefficients from G adj (τ ).Even though this correlation function has been studied in the past [49][50][51], its precise connection to quarkonium transport has remained unexplored, until now.It turns out that neither Eq. (10) nor Eq. ( 11) is valid for the quarkonium case.This is so because Eq. ( 10) is a result of the standard KMS relation, which, as we will show momentarily, is more complicated for the quarkonium correlator.Furthermore, Eq. ( 11) relies on the spectral function being odd in ω, which is crucially not true for the quarkonium correlator, as we will discuss in what follows.

KMS Relation and Non-odd Spectral Function
To explain the non-oddness of the spectral function for quarkonium transport, we follow Ref. [33] to use the usual proof of the KMS relation, plus the time-reversal operation and find which is the necessary KMS relation for proper thermalization of the internal degrees of freedom of the heavy quark pair (their relative motion and internal quantum numbers [52]).We then introduce the spectral function that governs quarkonium transport as which, by definition satisfies [g ++ adj ] > (ω) = (1 + n B (ω))ρ ++ adj (ω), with n B (ω) = (e βω − 1) −1 .We have kept the superscripts "++" in the label of this spectral function because we can also define which contains the same information, and satisfies Here comes the most important part: The spectral function (17) is not odd in ω.In the standard thermal field theory setup, we define ρ(ω) = g > (ω) − g < (ω) where g > (t) = ⟨ϕ(t)ϕ(0)⟩ and g < (t) = ⟨ϕ(0)ϕ(t)⟩, which are related via g > (ω) = g < (−ω) in frequency space by time translational invariance.This immediately leads to ρ(ω) = −ρ(−ω).However, the relation g > (ω) = g < (−ω) is not true for [g ++ adj ] > (ω) and [g −− adj ] > (ω) due to the path ordering of field operators and the additional Wilson line along the imaginary time in . Therefore, we do not know how ρ ++ adj (ω) transforms under ω → −ω a priori.
To see this more formally, one may also write the spectral function as a spectral decomposition in terms of the eigenvalues/eigenstates of H, denoted by {E n , |n⟩}, and those of [H1 − gA c 0 (0)T c adj ] ab , denoted by { Ẽn , |ñ a ⟩}, where a is interpreted as a component of the state, rather than a label.With these definitions, it follows that There is no reason why this expression would be odd under ω → −ω, because the energies E n and Ẽn can (and will) be different in general.Indeed, explicit perturbative calculations at NLO show that ρ ++ adj (ω) contains both ω-odd, which is the usual case, and ω-even parts (see Appendix B).The final result Eq. (3.66) shown in Ref. [33] is only for ω > 0, as mentioned there.We have performed a similar calculation for ω < 0 and found an ω-even part, which originates from the diagrams (5, 5r) of Ref. [33], or diagrams (j) of Refs.[34,45] where we have also added a factor of 2 since the definition of the spectral function shown in Eq. (3.66) of Ref. [33] differs from Eq. ( 17) by a factor of 2 (see Eq. (3.28) therein).
To demonstrate the importance of the ω-even part, we use it to recompute the difference between γ fund and γ adj at the order of α 2 where 2Nc .This difference was first calculated in Ref. [45].Some algebra and use of the definitions for [g ±± adj ] > leads to where we have used [g ±± adj ] > (ω) = (1+n B (ω))ρ ±± adj (ω) and used that they are translationally invariant in time.The piece proportional to θ(ω) is a pure vacuum contribution that vanishes in dimensional regularization.The second term inside the integral, however, gives a thermal contribution: which is exactly the difference given in Eq. ( 21).This settles a long-standing issue regarding the consistency of the gauge-invariant chromoelectric correlators in the adjoint and fundamental representation, and verifies explicitly that the spectral function relevant for quarkonium transport is qualitatively different from that for heavy quark diffusion.The above discrepancy ∆γ is explained precisely because ρ ++ adj (ω) is not odd in frequency.With these theoretical foundations in hand, we can now proceed to write down the formula analogous to Eq. ( 11), which will allow for the extraction of κ adj and γ adj from the evaluation of the Euclidean correlator G adj (τ ).

Extraction Formulas
Using the fact that G adj (τ ) is the analytic continuation of [g ++ adj ] > (t) to Euclidean signature, we can write However, in contrast to Eq. ( 11), the integrand may not be symmetrized with respect to ω because ρ ++ adj (ω) is neither even nor odd.We note that, as one might suspect from Eq. ( 15) and is apparent from Eq. ( 24), the analytic continuation holds provided that 0 < τ < β.This is precisely the range where we discuss the calculation of G adj in the next section.A direct calculation using Eqs.(11), (20), and (24) shows that where ζ(s, a) = ∞ k=0 (k + a) −s is the Hurwitz zeta function.
After extracting ρ ++ adj (ω) from the lattice QCD calculated G adj (τ ), which will be discussed in the next section, we can obtain κ adj and γ adj as where the expression we have written for κ adj makes it manifest that only the ω-odd part of ρ ++ adj (ω) contributes to it.(One can show this by using Eqs.( 1) and (4).)We note that γ adj may be substantially more difficult to extract than in the fundamental representation case.While the first term is indeed the same as in the fundamental case by virtue of , the fact that ρ ++ adj is not necessarily odd under ω → −ω means that the last term can contribute.Indeed, it does so in perturbation theory, as demonstrated by our calculation of ∆γ in Eq. (23).There is even an additional complication in that the 1 in 1 + 2n B of the second line will usually generate ultraviolet divergences that have to be regulated analytically (e.g., by dimensional regularization).Furthermore, the first term may also require regularization for the integration regions where τ ≈ 0, β.

IV. LATTICE QCD DETERMINATION OF G adj (τ ) AND RENORMALIZATION
In this section, we discuss how to perform a lattice QCD calculation of G adj and extract ρ ++ adj .We will first show a discretized version of G adj and then discuss how to renormalize the lattice QCD result when taking the continuum limit.Finally we will give a fitting ansatz to extract ρ ++ adj from the calculated G adj , which can then be plugged into Eq.( 26) to obtain the quarkonium transport coefficients.

A. Lattice Discretization
The main ingredient we require in order to construct a lattice formulation of the correlator that determines quarkonium transition rates is a discretized formula for the gauge field strength This discretization is different from the standard square plaquette.We chose this one because it makes the operator symmetric around the Wilson line direction, as shown in Fig. 1.One can then write an expression purely in terms of link variables for the correlator: where τ = an τ , and the products are ordered in such a way that the lower limit of the index labels corresponds to the operator that is most to the right in the product, and the upper limit to the one that is most to the left.A graphic representation of the correlator can be found in

B. Renormalization and Infrared Renormalon
The bare chromoelectric correlator G adj (τ ; a) can be evaluated by the lattice method explained above.For physical quantities, the lattice calculation result needs proper renormalization.Since the operator involves a Wilson line, it is expected that G adj (τ ; a) contains a linear divergence (which has not been explicitly checked and should be done so in the future via, e.g., a calculation in lattice perturbation theory), in addition to the usual logarithmic divergence.Therefore, we renormalize the bare correlator via where Z stands for the renormalization factor for the logarithmic divergence of the composite operator, with µ the renormalization scale and δm the mass renormalization associated with the self energy of the Wilson line.
It has been shown that this form of the renormalization factor for the nonlocal operator is consistent with the fact that when the nonlocal operator is expressed as a weighted sum of local lattice operators, they mix in the renormalization group flow [53].In this work, we will not address the potential mixing between similar correlators with different Wilson line paths connecting the two chromoelectric fields.
A NLO calculation of the real-time partner of G adj , i.e., [g ++ adj ] > has shown that [33] where we used Z ′ to distinguish the renormalization factor for [g ++ adj ] > from the Z for G adj .The "0" coefficient of the 1/ϵ term emphasizes that [g ++ adj ] > has no logarithmic divergence at NLO.The calculation was performed in the continuum by using dimensional regularization.The divergent term should be the same in the dimensionally regularized and lattice regularized perturbative calculations.Only the finite terms can be different.If we want to obtain the renormalized result in the MS scheme, the finite difference between the lattice scheme result and the MS result should still be accounted for.In the case of G fund , the difference is known at NLO [54].We leave the calculations of Z for the Euclidean G adj in both schemes to future studies.(As can be seen by comparing to Ref. [54], such calculations are research projects on their own.) Since the δm term is associated with the self energy of the Wilson line, one can use lattice perturbative calculations to determine it, but the uncertainties are expected to be large due to infrared renormalons.In particular, δm is expected to be of the form where m −1 is constant at leading order in lattice perturbation theory, but it has a residual dependence on a at higher orders via, e.g., aΛ QCD due to renormalization effects.On the other hand, m 0 is independent of the lattice spacing a, but it is scheme dependent as well.(Both m −1 and m 0 also depend on the other mass scales of the theory, if there are any.)The infrared renormalon ambiguity leads to an uncertainty in summing the perturbative series for m −1 , which is compensated by the same uncertainty in determining m 0 .The fact that both m −1 and m 0 are scheme dependent is reflected in the systematic uncertainty of fitting the a dependence from lattice calculations at small a, as shown in the recent study on renormalizing the quasi parton distribution function (quasi-PDF) [55].
Here we discuss a strategy to reduce the uncertainty caused by the infrared renormalons in determining the renormalization factor δm by using lattice QCD calculation results, which is motivated by the recent work on self renormalization of the quark quasi-PDF [55,56].The first step is to fit m −1 from the a dependence of Z(µ, a)G adj (τ ; a) when a is small for some τ .Different choices of τ are expected to give the same fitting result, as long as we maintain τ ≫ a to have negligible lattice artifacts).Due to the unknown nonperturbative dependence of m −1 on a, different parametrizations may be used in the fitting and they do not lead to the same result necessarily, which reflects the scheme dependence of m −1 .Then we define G R ′ adj (τ, µ) ≡ Z(µ, a)e m−1τ /a G adj (τ ; a), i.e., we only absorb the extracted a-dependent linear divergence and the logarithmic divergence into the renormalization factor and perform an operator production expansion (OPE) at small τ (i.e., β ≫ τ but we still require τ ≫ a) where O n denotes the local operators in the OPE and ⟨O n ⟩ R T (µ) represents their renormalized expectation values at the same temperature T .The expectation values of O n can be calculated by standard lattice QCD methods and renormalized perturbatively by calculating the corresponding logarithmic renormalization factors via lattice perturbative calculations, in the same way as it is done for the logarithmic renormalization factor Z for G adj .These local operators do not involve Wilson lines and thus do not have linear divergence, so it is expected that their renormalization is insensitive to the effects from infrared renormalons.The local OPE operators that may contribute include where e ρ is a unit vector along the spacetime direction ρ.The short-distance Wilson coefficients C n can be calculated in perturbation theory at the scale µ = 1/τ .The calculation of these coefficients is an active area of research [57,58].In practice, we can determine m 0 via Eq.( 32) by calculating the lattice renormalized G R ′ adj (τ, µ) and ) by including the renormalization factor associated with m 0 .As suggested in Ref. [56], to reduce the uncertainty caused by the infrared renormalons, one resums the leading infrared renormalons in C n by regulating the renormalon poles in the Borel space and applying the inverse Borel transformation.As shown therein, this strategy removes a large uncertainty in the determination of the quark PDF.We expect a similar uncertainty reduction to happen for the determination of G R adj by using this strategy.After determining the renormalized G R adj in the lattice regularization, we can convert it into the MS scheme if we know the difference between the perturbative results of the logarithmic divergence in these two schemes.As part of the conversion process, one has to take care of the fact that in dimensional regularization with d = 4 − ϵ and ϵ → 0, the linear divergence is absent.Any residual finite terms from this linear divergence are accounted for through m 0 in the OPE matching.

C. Fitting Ansatz for ρ ++ adj
Once we obtain the renormalized G R adj , we can use Eq. ( 24) to fit the spectral function ρ ++ adj .Since we only have a limited number of data points in τ , we need a fitting ansatz.One ansatz that has been used in the lattice studies of the heavy quark diffusion coefficient is of the form [46] where ρ IR and ρ UV represent ansatzes in the small and large ω regions, respectively.We will construct ansatzes motivated from perturbative studies.
The ω-even part of the large frequency behavior of ρ ++ adj (ω) is determined by Eq. (20).The remaining (ωodd) terms can be read off directly from [33], where ρ ++ adj (ω) was calculated at ω > 0. Explicitly, where N f is the number of light (massless) quark flavors in the theory.It was shown in Ref. [34] that up to O(g 4 ), the leading temperature-dependent contributions (which are the same for ρ ++ adj and ρ fund , cf. [33]) at large frequency go as T 4 /ω, which are omitted in Eq. ( 35) since they are subleading.
On the infrared side, one needs to use the hard thermal loop effective theory to capture the behavior of correlation functions when |ω| ≲ gT ∝ m D , where m D is the socalled Debye mass of the QGP, given (perturbatively) by , which quantifies color-electric screening in a thermal plasma.To see the difference between the ρ fund and ρ adj in the small ω region, one needs to consider the same type of diagrams that led to the difference shown in Eq. ( 20), which has a prefactor of g 4 , meaning that the dominant corrections in the regime |ω| ≲ m D will be of order g 4 m 2 D |ω| ∝ g 6 T 2 |ω|.This means that we cannot make quantitative statements by considering only the 1-loop diagram that leads to Eq. ( 20) (replacing the propagators with their HTL-resummed counterparts), as we can get competing effects from 2loop diagrams in QCD, which contribute at order g 6 .In practice, one would also need to calculate these 2-loop diagrams to be able to match the HTL result to full QCD.We will leave such calculations to future studies.Here we only list the leading contribution in the infrared regime, which can be written in terms of the well-known heavy quark diffusion coefficient κ fund at NLO: where κ fund is given by [36,59]: with C ≈ 2.3302, as given in Ref. [59].The fact that the low-frequency limit of the adjoint and fundamental correlators do not differ up to this order had already been noticed in Ref. [59].Motivated by the above perturbative analyzes, we suggest to use Eq. ( 35) as ρ UV in the fitting ansatz (34) and use κ adj ω + c|ω| to parametrize ρ IR with c some constant that does not contribute to κ adj .The appearance of the c|ω| term in ρ IR is a crucial difference from the case of the heavy quark diffusion coefficient and is motivated by perturbative calculations shown in Section III C 1.The fitting of ρ ++ adj will not only provide the quarkonium transport coefficient κ adj , but also the frequency dependence of ρ ++ adj , which is important to evaluate γ adj , as well as the frequency-dependent correlators g ±± adj (ω) that determine the quarkonium dissociation and recombination rates.

V. CONCLUSIONS
In this paper, we explained how to determine the real time quarkonium transport properties from a Euclidean chromoelectric field correlator.This determination requires to reconstruct a spectral function in a way that is different from more intensively studied spectral function reconstruction problems, such as the one required for the extraction of the heavy quark diffusion coefficient.The key results are shown in Eq. ( 26).We then discussed the lattice determination of the Euclidean correlator, and in particular, a method to reduce the uncertainty caused by infrared renormalons in obtaining the renormalization factor for the linear divergence of the correlator.This method is quite involved and several perturbative calculations needed to implement the method are left to future studies, such as the lattice-regularized perturbative calculation of the logarithmic renormalization factor Z in Eq. ( 29) and the Borel-resummed calculation of the Wilson coefficients in the OPE (32).Our work paves a way towards a nonperturbative determination of the quarkonium transport properties in the QCD hot medium, which generalizes the use of a weakly interacting gas of quarks and gluons as a microscopic model of the QGP in Boltzmann (rate) equations [60][61][62][63] for quarkonium to the strongly coupled case.This not only deepens our understanding of the QGP and quarkonium production in heavy ion collisions, but may also provide insights for studies of exotic heavy flavor production [64][65][66] and dark matter bound state formation in the early universe [33,[67][68][69].
After publication of this work, we noticed that although the difference ∆ρ between spectral functions we found in this work in Eq. ( 20) reproduces the difference between the transport coefficients γ adj and γ fund present in [45], it does not reproduce the difference between the coefficients of the π 2 terms in the result of the calculation for ρ ++ adj in [33] and the result of the Euclidean QCD calculation of ρ fund in [34].In the original version of this work, we made a transcription error that led us to find no tension.It is clear from earlier calculations by others [49] and ourselves [33] that the ω > 0 part of the spectral function ρ ++ adj at T = 0 is as given in Eq. (35).It is also clear that the only ω-even contributions to ρ ++ adj are given by Eq. ( 20), as ρ fund is odd in ω and Eq. ( 20) unambiguously represents the subtraction of the ω-odd contribution from diagrams ( 5) and (5r) of [33].Therefore, the ω < 0 part of ρ ++ adj at T = 0 can be cross-checked with the above two calculations and is as given in Eq. (35).There is no thermal contribution to the difference at 1loop because this difference can be written in terms of the imaginary part of free retarded correlators, which do not depend on T .To unambiguously establish the origin of the discrepancy between the π 2 terms in ρ ++ adj and ρ fund , and as an additional cross-check, we have calculated the difference between G adj and G fund in the imaginary time formalism and added it as Appendix C to this work.We find that the difference between ρ ++ adj and ρ fund is as given in Eq. ( 20), and therefore, that the term proportional to π 2 in ρ fund should be −π 2 /6 instead of −2π 2 /3.That is to say, we find that in [34], −8π 2 /3 in Eq. (4.2) should be −2π 2 /3.This new calculation will appear as part of the erratum in the journal version.The explicit calculation of this integral is equivalent to the one presented in the Supplemental Material of Ref. [35].The final result is as claimed in the main text.
It is noteworthy that the difference between the spectral functions, as given in Eq. (B2) may also be used in conjunction with HTL-resummed propagators to explore the value of the difference (a modification to the gluon 3-vertex is also necessary, according to the HTL effective theory Feynman rules.They can be found in Ref. [70].).However, as discussed in the main text, a full fixed-order calculation at O(g 6 ), which is the leading contribution to the difference in the small frequency domain, also requires considering 2-loop diagrams, which we will not pursue here.
It then follows that the difference between spectral functions in the limit D → 4 is purely determined by the divergent contribution to ∆ G(k n ).In the limit, (4 − D)Γ(4 − D) → 1, and we may set D = 4 elsewhere.The integral over Feynman parameters gives a simple result just as we obtained via our real time calculation.

C.1. Resolving the tension with previous results
We now discuss the calculation in [34] of the terms proportional to π 2 , and how they arrived at a different result.In short, the issue is that the IR regulators employed in [34] fundamentally alter the analytic structure of the integral to be calculated.
In [34], the integral structure that generates the π 2 terms is given by their Eq.(A.42) 2. Graphic representation of the pole structure of Eq. (C7) in the complex kn plane at fixed k.Poles are represented with blue crosses, the branch cut between −ik and +ik is represented by a red zig-zag line, and the branch cuts above and below ±ik are represented by wavy green lines.The latter branch cut is induced by the presence of the Wilson line.Left: The pole structure without a regulator.Right: The pole structure with the regulator used in [34].Without the branch cuts induced by the Wilson line, this regulator is not problematic because the branch cut denoted by a red zig-zag line does not intersect the poles.However, with the branch cuts induced by the Wilson line, moving the poles in this manner qualitatively alters the pole structure, because the contributions from the regions where ik < ±Im(kn) < i √ k 2 + λ 2 , Re(kn) ≈ 0 will contribute with an opposite sign to the unregulated version.

Fig. 1 .
The average ⟨•⟩ E represents the expectation value under the measure defined by the Euclidean lattice path integral, i.e., ⟨O⟩ E = 1 Z E DU exp(−S E [U ])O[U ] where Z E = DU exp(−S E [U ]).

FIG. 1 .
FIG. 1. Lattice discretization of the chromoelectric field correlator.The electric field insertions are constructed by taking the difference between the products of gauge links over the blue and red contours at the ends of the light blue contours, which represents an adjoint Wilson line.In this setup, the adjoint Wilson line is equivalent to two antiparallel fundamental Wilson lines.