Chromoelectric field correlator for quarkonium transport in the strongly coupled $\mathcal{N}=4$ Yang-Mills plasma from AdS/CFT

Previous studies have shown that a gauge-invariant correlation function of two chromoelectric fields connected by a straight timelike adjoint Wilson line encodes crucial information about quark-gluon plasma (QGP) that determines the dynamics of small-sized quarkonium in the medium. Motivated by the successes of holographic calculations to describe strongly coupled QGP, we calculate the analog gauge-invariant correlation function in strongly coupled $\mathcal{N}=4$ supersymmetric Yang-Mills theory at finite temperature by using the AdS/CFT correspondence. Our results indicate that the transition processes between bound and unbound quarkonium states are suppressed in strongly coupled plasmas, and moreover, the leading contributions to these transition processes vanish in both the quantum Brownian motion and quantum optical limits of open quantum system approaches to quarkonia.


Introduction
Heavy quarkonium is a bound state containing a heavy quark-antiquark pair (QQ). Due to the large heavy quark mass, the spectra of the ground and the first few excited quarkonium states can be studied by solving a Schrödinger equation with a potential model such as the Cornell potential [1]. When a quarkonium state is placed inside a hot and/or dense nuclear environment, such as the quark-gluon plasma (QGP) produced in heavy ion collisions, the attractive potential can be significantly screened. As a result, at high temperatures a heavy quark-antiquark pair can no longer form a bound state, i.e., quarkonium melts [2,3]. This static screening picture has motivated using quarkonium as a signature of the QGP formation in heavy ion collisions and, more generally, as a probe of the QGP. However, dynamical processes such as dissociation [4][5][6][7][8][9] and recombination [10][11][12] complicate the description of quarkonium inside the QGP. At the lowest order in a weak coupling picture, dissociation occurs when a quarkonium state absorbs a gluon from the medium and the bound state is excited to the continuum [8], while recombination happens when an unbound QQ pair radiates out a gluon that takes away enough energy so that the unbound pair forms a bound state [13]. Many studies have been devoted to calculate these dynamical processes in a weak coupling picture, which can be dated back to the early work by Bhanot and Peskin [14,15]. However, it is known from current heavy ion collision experiments that the QGP created in such collisions is actually a strongly coupled fluid, which means that calculations derived from a weakly coupled description are not always reliable for practical phenomenology. More specifically, the issue of a weak coupling picture is that the QGP medium relevant for quarkonium dissociation and recombination is far away from being a gas of weakly interacting light quarks and gluons.
The question then becomes how quarkonium dissociation and recombination can be formulated for a strongly coupled medium, where perturbative techniques may not apply. Thanks to recent developments using the open quantum system framework and nonrelativistic effective field theories of QCD, both quantum and classical transport equations to describe quarkonium in-medium dynamics have been derived  (see recent reviews [40][41][42][43]). These developments allow one to factorize out the nonperturbative soft physics of the QGP from the rest of the heavy quark physics [44]. More specifically, the nonperturbative soft physics that originates from the strongly coupled QGP is encoded in terms of a gauge-invariant correlator of chromoelectric fields dressed with Wilson lines [44] (see also Ref. [23,27] for a different perspective). So far, this chromoelectric field correlator is only known up to next-to-leading (NLO) order in perturbation theory at finite temperature [45], which already shows a difference from a similar but different chromoelectric field correlator that characterizes single heavy quark diffusion [46,47]. In spite of the fact that this correlation function has been previously considered in vacuum and Euclidean QCD (see, for instance, Refs. [48][49][50][51][52][53][54]), not much is known about the chromoelectric field correlator for quarkonium in-medium dynamics beyond NLO. With the proper setup, one may be able to calculate the zero frequency limit of the correlator by using Euclidean lattice QCD methods and analytically continuing the result to real time. Studying the correlator at finite frequency directly in real time is intractable numerically due to the notorious sign problem.
In this paper, we calculate this chromoelectric field correlator in strongly coupled N = 4 Yang-Mills theory at both zero and finite frequency using the AdS/CFT correspondence [55,56]. The general strategy of our studies parallels the approaches used in studies of jet quenching and heavy quark diffusion, where effective field theory was first applied to describe the strongly coupled physics inside the QGP in terms of gauge-invariant objects that require a nonperturbative determination, which were then subsequently calculated via the AdS/CFT holographic duality. In the case of jet quenching, Soft-Collinear Effective Theory [57][58][59][60][61] had been used to formulate the jet quenching parameter as a Wilson loop consisting of two light-like Wilson lines, and was later calculated using the AdS/CFT correspondence in Refs. [62][63][64]. In the case of heavy quark diffusion, heavy quark effective theory was applied to define the heavy quark diffusion coefficient as a chromoelectric field correlator in Ref. [65], which was then calculated in the same work via the holographic principle (see Ref. [66] for a different way of studying heavy quark diffusion in the AdS/CFT approach). In this work, we will use the AdS/CFT technique to calculate the chromoelectric field correlator relevant for quarkonium transport, which has been defined by using potential nonrelativistic QCD (pNRQCD) [67,68]. This correlator is given by [23,27,44,45,47] T E a i (t)W ab where E a i is the chromoelectric field, W ab [t,0] is an adjoint representation Wilson line andT denotes the time ordering operator. Our approach will be to first rewrite the chromoelectric field correlator as a path variation of a Wilson loop, which defines a contour on the boundary of an AdS black hole spacetime. Then, using the holographic correspondence, we will calculate the expectation value of the Wilson loop by finding the extremal surface in the bulk of the AdS spacetime that hangs from the contour defined by the Wilson loop. Finally, we will obtain the expectation value of the chromoelectric field correlator by taking the path variation, which amounts to solving linear equations for fluctuations that propagate on the extremal surface.
The result presented here is important for quarkonium phenomenology, because it provides the first nonperturbative picture of the in-medium quarks and gluons that are relevant for quarkonium dynamics in the QGP and goes beyond the assumption of a weakly interacting gas. Crucially, the nonperturbative distributions of in-medium quarks and gluons can be process-dependent, as in the case of deep-inelastic scattering (inclusive versus semi-inclusive, polarized versus unpolarized, and so on). In general, it is not expected that the in-medium quarks and gluons relevant for jet quenching would have the same distribution as those affecting heavy quarks.
The paper is organized as follows: We will first review the transport equations for quarkonium in-medium dynamics and the relevant chromoelectric field correlator in Section 2. The setup of the AdS/CFT calculation will be given in Section 3, followed by details of the calculation in sections 4 and 5 for two well-motivated setups. We will discuss the AdS/CFT result with a focus on its implications for quarkonium transport in the QGP in Section 6. Finally, we will draw our conclusions in Section 7.
2 Quarkonium transport at high and low temperatures Three energy scales are used to describe heavy quarkonium in vacuum: the heavy quark mass M , the inverse of quarkonium size 1 r and the binding energy |E b |. Nonrelativistically, these three scales form a hierarchy M ≫ 1 r ≫ |E b |. Depending on where the plasma temperature T fits into the hierarchy, we may have different descriptions of quarkonium in-medium dynamics. We will always consider the case of M ≫ T , since the highest temperature achieved in current heavy ion collision experiments is on the order of 500 MeV, which is smaller than the charm quark mass M c ≈ 1.3 GeV and much smaller than the bottom quark mass M b ≈ 4.2 GeV. Furthermore, throughout this paper, we will focus on heavy quarks and quarkonia at low transverse momenta.
When the plasma temperature is very high M ≫ T ≫ 1 r ≫ |E b |, which is the case in the early stage of heavy ion collisions, the interaction between a heavy quark-antiquark (QQ) pair is significantly screened. As a result, the in-medium dynamics of a QQ pair can be described in terms of two independent heavy quarks that diffuse and dissipate in the plasma. This dynamics can be approximately described by a Langevin equation with drag and diffusion and the heavy quark diffusion coefficient has been calculated nonperturbatively via lattice methods [69][70][71][72]. One can improve the Langevin equation by adding an attractive potential effect that is screened when the QQ pair is far away in space. This potential effect can only last for a time scale given by the imaginary potential obtained from lattice studies [73][74][75][76][77].
After the formation of the QGP in heavy ion collisions, the plasma expands quickly and cools down. When the temperature drops to the region where T ∼ 1 r , the interaction between a QQ pair can no longer be neglected and must be included in the Langevin description. When the temperature further drops, M ≫ 1 r ≫ T ≫ |E b |, a different description comes into play, in which the color correlation between the QQ pair becomes important. The in-medium dynamics of a QQ pair can be described by a Lindblad equation in the quantum Brownian motion limit [23,27] (the subscript "S" stands for "system" in the open quantum system framework, as contrary to the environment, i.e., the QGP) where the Hamiltonian is given by , (2.2) and the density matrix is assumed to be block diagonal in the color singlet and octet basis ρ S (t) = ρ (s) . (2. 3) The Lindblad operators are given by where i = x, y, z. The transport coefficients κ adj and γ adj are defined in terms of chromoelectric field correlators where T F is defined by Tr(T a F T b F ) = T F δ ab with T a F being the generator in the fundamental representation and W ab (t, 0) denotes a time-like Wilson line in the adjoint representation from time 0 to t: W(x, y) = P exp ig x y dz µ A a µ (z)T a A , (2.9) in which x and y are Minkowski position 4-vectors connected by a straight path. The expectation value of an operator O is defined as ⟨O⟩ T ≡ Tr(Oe −βH E )/ Tr(e −βH E ) where β = 1/T is the inverse of the plasma temperature and H E denotes the Hamiltonian of light quarks and gluons in the QGP. As can be seen, in this high temperature limit where the quantum Brownian motion limit of the open quantum system framework is valid, it is the zero frequency limit of the correlator shown in Eq. (1.1) that is relevant for quarkonium in-medium dynamics.
As the QGP temperature continues dropping and finally becomes of the same order as the binding energy, M ≫ 1 r ≫ T ∼ |E b |, we need another description that is based on a classical Boltzmann equation which can be derived by using the open quantum system framework in the quantum optical limit, pNRQCD, the Wigner transform and semiclassical gradient expansion [44] (a subtlety of using the quantum optical limit can be resolved by working in the semiclassical limit, as explained in Ref. [42]). If we further integrate over the momentum distribution of the phase space distribution, we will arrive at a rate equation for the density of a quarkonium state n b with the quantum number b, 10) where Γ is the dissociation rate and F denotes the contribution of quarkonium formation (in-medium recombination). They are given by x, p cm , x rel = 0, p rel ) , (2.12) where ⟨ψ b |r|Ψ p rel ⟩ is the dipole transition amplitude between a bound quarkonium state wavefunction ψ b and an unbound QQ wavefunction Ψ p rel that is a scattering wave with momentum p rel . The two-particle phase space distribution f QQ (t, x, p cm , x rel = 0, p rel ) is for an unbound QQ pair with the center-of-mass (cm) position x cm = x, cm momentum p cm , relative position x rel = 0 and relative momentum p rel . The relative position is fixed to be 0 which is a result of a gradient expansion used in taking the semiclassical limit.
The QQ phase space distribution does not factorize into the product of two single particle distributions f QQ (t, x, p cm , x rel = 0, p rel ) ̸ = f Q (t, x, p Q )fQ(t, x, pQ) , (2.13) which means that the formation term F in the rate equation can account for both correlated and uncorrelated recombination [78]. The Boltzmann and rate equations have been extensively used in phenomenological studies of quarkonium and exotics production in heavy ion collisions [79][80][81][82][83][84].
In the rate equation, we follow the notation of Ref. [45] and write the physical (Wightman) correlation functions that govern quarkonium transport as Tr H E a i (t)W ac (t, +∞)W cb (+∞, 0)E b i (0)e −βH (2.14) [g ++ E ] < (t) = 1 Z Tr H W cb (+∞, 0)E b i (0)E a i (t)W ad (t, +∞)e −βH W dc (+∞ − iβ, +∞) (2.15) where H is the environment (QGP) Hamiltonian, Tr H denotes a trace over states in the Hilbert space of the theory, and Z = Tr H e −βH is the partition function of the QGP. As opposed to Ref. [45], here we have written the thermal averages explicitly, with the (Euclidean) adjoint Wilson line at Re{t} = −∞ in the definition of [g −− E ] > (t) accounting for the fact that in the corresponding physical situation, the QGP and the point adjoint color charge have thermalized together.
In frequency space, we define the correlators as 18) in terms of which the KMS relations are given by 20) and are related by the combined action of parity and time-reversal where we have assumed the QGP is symmetric under these discrete symmetries. A short proof of the KMS relations is given in Appendix A, where we give a derivation that is more clearly indicative of the physics at work than in Ref. [45]. From this proof it is also clear that the normalization factor Z of the KMS conjugates [g ±± ′ E ] ≶ (ω) must be the same. Perturbative calculations [45] verify that the QGP partition function adequately serves this purpose.
Following [45], one may then define spectral functions that respect the KMS relations: Contrary to typical spectral functions, the ones we have just introduced are not guaranteed to have a definite parity under ω → −ω. Rather, because of parity and time-reversal, they are related to each other via: so it is still true that complete knowledge of one spectral function fully determines all correlation functions introduced above. Conversely, complete knowledge of one of the physical (Wightman) correlation functions also fully determines all the rest of the correlation functions.
In this paper, we will use the AdS/CFT correspondence to evaluate the time-ordered version of the correlation function: which can be written as (t) in hand we also have access to the anti-time-ordered correlator through complex conjugation (denoted by a star * ) Then, we can also express the Wightman function [g ++ E ] > (t) in terms of its (anti)timeordered counterparts: This means that once we obtain the time-ordered correlator from the AdS/CFT correspondence, we can evaluate all the other correlation functions, in particular the physical (Wightman) correlators that enter the Boltzmann and rate equations. We can use all of the above to write Eq. (2.28) in frequency space: where P denotes the Cauchy principal value. The inverse map is given by 1 (ω), the spectral functions follow immediately. As such, all we need is to calculate the time-ordered correlator. In the next section we describe the theoretical foundations we will employ to evaluate it.

Field strength correlators from Wilson loops
In this section we will describe the setup of the calculation of the non-Abelian electric field correlator we wish to obtain. We will start in Subsection 3.1 by describing how such a correlation function can be obtained by taking variations of a Wilson loop in a purely fieldtheoretic setup. Then, in Subsection 3.2 we will proceed to describe how we can evaluate Wilson loops at strong coupling in N = 4 supersymmetric Yang-Mills (SYM) theory using the AdS/CFT correspondence. Here we will discuss the role of the additional parameter n ∈ S 5 when constructing the supersymmetric Wilson loop that preserves the features of an adjoint Wilson line. The subsequent Subsection 3.3 discusses how to take variations of a Wilson loop (i.e., how to introduce field strength insertions) along the contour that defines it from the point of view of the dual gravitational theory. We will also establish the prescriptions necessary to fully define the correlation function for quarkonium in-medium dynamics, and discuss the differences with the correlation function that determines the heavy quark diffusion coefficient.
We will motivate two different holographic configurations that determine different correlation functions, and describe the setup for the calculation of each object in separate sections 4 and 5. As we will shortly discuss, the configuration we present in Section 5 possesses features that make it a better proxy for the correlation function we seek to calculate (3.12) than the one in Section 4. Nonetheless, to be thorough, we will conclude which correlator is a better N = 4 SYM proxy for the description of two non-Abelian electric fields connected by an adjoint Wilson line only after obtaining explicit results for both of them. 1 To prove that they are the inverse of each other, one needs to use a particular representation of the Dirac delta: which may be verified by direct action of this distribution on functions whose arguments are the parameter a.

Wilson loops in gauge theory and their variations
As we pointed out in the previous section, the transport properties of quarkonia are governed by the correlation functions of chromoelectric fields dressed by Wilson lines, calculated inside a medium. Properties of the medium are encoded in terms of expectation values, which can be determined by a thermal density matrix Z −1 exp(−βH). However, for the purposes of this subsection we will not need to make reference to the nature of the medium. Rather, we will only discuss how to construct the expectation value we are interested in, by starting from another class of observables that has a well-defined prescription for evaluation at strong coupling using the gauge-gravity duality, which we will discuss explicitly in Subsection 3.2.
Concretely, it is possible to study the correlation functions of gauge theory field strengths F µν dressed by Wilson lines starting from another class of gauge-invariant operators, namely, Wilson loops W [C]. They are defined as where A a µ is the SU(N c ) non-Abelian gauge potential, T a denotes the generator matrix of the group, g is the coupling constant, and P denotes path ordering in the product of group elements A µ = A a µ T a . The study of these operators is a cornerstone of much of our understanding of heavy or highly-energetic quarks, and in particular for their dynamics inside a thermal medium. The heavy quark-antiquark interaction potential [85][86][87], the jet quenching parameter [62,88], and the heavy quark diffusion inside a thermal medium [65,89,90] have all been formulated and studied through Wilson loops. Crucially, all of them admit a holographic description in N = 4 supersymmetric Yang-Mills theory.
As it turns out, one can connect the expectation values of Wilson loops and that of field strengths dressed by Wilson lines by considering functional variations of the path on which the Wilson loop is defined. Concretely, we consider a Wilson loop defined by a path C, and let γ µ (s) be a parametrization of the path, with s ∈ [0, 1], and γ µ (0) = γ µ (1) for a closed path. Then we consider a deformation of the path, which we denote by C f and parametrize by γ µ f (s) = γ µ (s) + f µ (s). It is then a textbook exercise [91] to show that where U [s ′ ,s] denotes a Wilson line from γ(s) to γ(s ′ ) in the same representation as W [C f ] andγ ρ (s) = dγ ρ (s)/ds. We will abuse the notation a bit to use U [γ(s ′ ),γ(s)] and U [s ′ ,s] interchangeably. The order of the operators in this expression, and of those in the present discussion before Section 3.2, only refers to the SU(N c ) matrix product. Operator ordering in the sense of the order in which they act on a state in the Hilbert space of the theory has yet to be specified (this will be done below, and further developed in Appendix B). This property provides us with a tool to generate as many field strength insertions as we want along the path of the Wilson loop. By acting on W [C] with one more derivative and assuming s 2 > s 1 , We see explicitly that we can obtain correlation functions of field strength operators F µν dressed with Wilson lines by taking derivatives with respect to the path on which the Wilson loop (3.1) is defined. 2 It is of course possible to continue beyond two field strength insertions, but for our present purposes it is sufficient to evaluate the two-point deformations. Specifically, our correlator of interest can be obtained by taking C to be a closed loop parametrized by γ µ (s), where s ∈ [−T /2, 3T /2], as with t µ = (1, 0, 0, 0) being a unit vector in the positive time direction. This describes a timelike loop that backtracks upon itself after reaching a maximal value for the time coordinate t = T /2. We note that in our setup we must have W [C f =0 ] = 1, even for time-ordered operators (see Appendix B for details). Then, taking variations of the path with respect to perturbations in a spatial direction f i (s) leads to where we have assumed that s 2 > s 1 . For our purposes, we will take all quantum mechanical operators to be time-ordered, which is exactly what one obtains from the path integral formulation of QFT. Some comments on the operator ordering can be found in Appendix B, where we discuss similarities and distinctions with other types of ordering, as well as explain why the time-ordered correlator describes quarkonium dynamics. 2 The only subtlety in this expression is that if s1 = s2, then there is another term on the right hand side, due to the action of δ/δf on theγ present in Eq. (3.2). This term is naturally proportional to derivatives of the delta function δ(s2 − s1), and can be easily isolated from the rest of the correlator by looking at positions s2 > s1 and continuously extending the result to s2 = s1 (whenever possible). If one looks at the Fourier transform of the left hand side of Eq. (3.3) with respect to s2 − s1, as we will find most natural to do later on, then there will be a contribution from points with s2 = s1 in the form of a polynomial of positive powers of their Fourier conjugate variable, which we will have to subtract to obtain the correlation function of interest.
We can further simplify this expression by restricting ourselves to contour deformations of the form where h µ (t) is the independent function with respect to which we will take a variation. Going through the same arguments given as above, one finds The factor of 4 is simply a consequence of doubling the size of the deformation to the contour. Intuitively, the reason why the correlator is determined by the antisymmetric deformation is that the symmetric contour deformation, i.e., f µ of the form with g µ (t) as an independent function, gives a vanishing variation which is a consequence of the more fundamental statement that W [C f ] = 1 for an arbitrary path change given by g µ (t) that can be even non-infinitesimal. The reason behind this is that two anti-parallel fundamental Wilson lines U, V that make up the Wilson loop as W = 1 Nc Tr color [U V ] will still satisfy V = U −1 , even when the path over which they are laid is deformed non-infinitesimally by f . On the other hand, the antisymmetric deformation described by h µ (t) gives a nontrivial result, precisely because in the language we just introduced, we have V ̸ = U −1 .
With all of this, we can state that the correlation function we want to calculate is given in terms of path variations of a Wilson loop through the following expression: where ⟨· · · ⟩ denotes the expectation value. Having formulated a way to obtain the desired correlation function in terms of operations upon Wilson loops, we now turn to the computational tool that we will use to evaluate this correlation function in N = 4 SYM theory.

Wilson loops in AdS/CFT
The AdS/CFT correspondence has proven to be an invaluable tool to gain insight into the strongly coupled regime of non-Abelian gauge theories, by casting a potentially intractable non-perturbative quantum mechanical problem for a conformal field theory (CFT) in terms of a purely classical problem in a concrete gravitational setup of higher dimensionality. In what follows, we will use the real-time formulation of the duality [56] between N = 4 supersymmetric Yang-Mills theory in the Minkowski four dimensional spacetime (Mink 4 ) in the large N c limit and type IIB string theory in an (asymptotically) AdS 5 ×S 5 spacetime. The asymptotic boundary of AdS 5 is identified as the Mink 4 on which the supersymmetric Yang-Mills theory lives. At strong coupling in the SYM theory, the dual description on the string theory side reduces to classical string dynamics in a curved spacetime.
The task now is to describe the expectation value of Wilson loops using the duality. This was first done by Maldacena in Ref. [86]. However, due to the supersymmetric nature of N = 4 SYM, the CFT object that has a simple gravitational dual description in AdS 5 is the locally 1/2 BPS 3 Wilson loop whereẋ µ (s) = dx µ (s)/ds and ⃗ ϕ = (ϕ 1 , . . . , ϕ 6 ) are the six Lorentz scalar fields in the adjoint representation of SU(N c ) intrinsic to N = 4 SYM. These scalars enter the Wilson loop coupled to a directionn(s) ∈ S 5 that specifies the direction along which the string (to be introduced in the next paragraph) "pulls" the heavy quark (the string lives in a 10 dimensional space AdS 5 × S 5 and has a string tension). It can be thought of as an additional property that the heavy quark carries as it propagates through Mink 4 , and in general, it depends on the coordinate s along the path C.
Specifically, the AdS/CFT duality gives an explicit prescription to evaluate the expectation value of these generalized Wilson loops. It is given by where is the Nambu-Goto action of a string configuration Σ described by X µ (τ, σ) ∈ AdS 5 × S 5 , with µ ∈ {0, 1, . . . , 9}. Σ(C;n) is the surface (also referred to as a "worldsheet") that extremizes the Nambu-Goto action S NG with Dirichlet boundary conditions given by C andn at the asymptotic boundary of AdS 5 . It is in this sense thatn defines the direction along which the string "pulls" the heavy quark. 4 The subtraction by S 0 [C;n] is necessary to regularize the result by subtracting the energy associated to the mass of the heavy quark propagating along C, and is also useful because by subtracting the rest mass of the heavy quark it isolates the "energy" associated to the interactions described by the Wilson loop at hand. (In QCD, the subtraction also contains the mass renormalization caused by the heavy quark self interaction that is part of the physics contained in the Wilson loop. However, in N = 4 SYM the self interaction diagrams are ultraviolet (UV) finite, and thus the subtraction only contains the bare heavy quark mass [92].) Only after this subtraction has taken place, may one have an action with a finite value that allows for a comparison of the "energy" of different string configurations determined by S NG . In the case where one may equate the expectation value of a Wilson loop to that of a time evolution operator over a time period T for a fixed set of boundary conditions C andn, it is the lowest energy E[Σ(C;n)] = (S 0 [C;n] − S NG [Σ(C;n)]) /T configuration that determines the expectation value mentioned above.
We also need to specify the background metric g µν that describes the dual AdS 5 × S 5 spacetime. Because of our interest in thermal physics, we use the metric for the dual description of a field theory at finite temperature T . In terms of Poincaré coordinates, it is given by [56] where f = 1 − (πT z) 4 , R is the curvature radius of the AdS metric, z is the radial AdS coordinate, with the asymptotic boundary at z = 0 and the black hole horizon at z = (πT ) −1 . The coordinates t and x describe Minkowski spacetime at z = 0, and the S 5 coordinates describe a sphere of radius R. This is the Schwarzschild-AdS metric that is dual to N = 4 SYM at finite temperature. The duality also prescribes R 2 /α ′ = √ λ = g 2 YM N c , and this single-handedly determines how the coupling constant appears in the results for the strongly coupled limit. In the language of Refs. [93,94], this corresponds to a single copy of a Lorentzian manifold in the gravity side of the duality. Therefore, all calculations done in this setup will give time-ordered quantities in terms of the action of operators on the Hilbert space of the quantum theory. If we needed to consider more complicated operator orderings we would be forced to introduce a larger manifold on the gravity side, containing more than one copy of AdS 5 × S 5 . This is manifest in the holographic descriptions of heavy quark diffusion and jet quenching [64,65], and is discussed at length in Ref. [93]. 4 To picture this, it is helpful to note that, asymptotically as z → 0, the part of the metric in Eq. (3.16) involving z andn is proportional to dz 2 + z 2 dΩ 2 5 , and as such, (z,n) may be thought of as a 6-component vector along the directionn with length z.

The role ofn
The classic results for Wilson loops in the strong coupling limit that are obtained from the gauge/gravity duality feature a constant value ofn throughout the quark trajectory [62,65,86]. This is a natural choice because, following the discussion of Ref. [86], the value ofn is determined by the vacuum expectation value of a Higgs field that has undergone spontaneous symmetry breaking. Given that this is tantamount to a choice of vacuum, the low-energy physics will not modify the value ofn throughout the trajectory of a heavy W-boson (which is later referred to as "quark," because the Wilson loop also describes the evolution of a heavy quark.). However, inspecting the expression for the 1/2 BPS Wilson loop (3.13), one realizes that keeping a constant value ofn =n 0 throughout both sides of the contour shown in Eq. (3.4) violates our expectation for conventional Wilson lines that a closed Wilson loop consisting of two overlapping anti-parallel Wilson lines has W [C] = 1. Indeed, the famous result for the heavy quark interaction potential [86], determined by evaluating a rectangular Wilson loop C L of temporal extent T and spatial size L (after subtracting the heavy quark mass contribution S 0 [C L ;n =n 0 ], and assuming T ≫ L) gives which does not satisfy lim L→0 W [C L ] = 1. The reason behind this is that the contributions from the scalar fields to the locally 1/2 BPS Wilson loop do not cancel ifn is held constant [86,92,95]. In other words, this way of approaching a Wilson loop made of coincident anti-parallel fundamental Wilson lines does not have the "zig-zag" symmetry [96,97] that the standard Wilson loop (3.1) respects. 5 This is an unavoidable obstacle if one attempts to interpret the variations of the Wilson loop that leads to the heavy quark potential as the dual object of a correlation function of two chromoelectric fields connected by an adjoint Wilson line. Concretely, the standard gauge theory identity for an adjoint representation Wilson line in terms of fundamental representation Wilson lines, i.e., Eq. (3.7), is not satisfied if one tries to build such an object by taking variations of a loop at constantn, because the two anti-parallel fundamental Wilson lines that enter the loop at constantn (U [t 2 ,t 1 ] and U [t 1 ,t 2 ] ) are not each other's inverse. This had already been noted and discussed previously in Refs. [86,92,95].
At this point, the appropriate Wilson loop to use in N = 4 SYM is apparent: We have to construct a locally 1/2 BPS Wilson loop that has two timelike links U, V of temporal extent T that satisfy V = U −1 , such that W BPS = 1 Nc Tr [U V ] = 1. This is realized when the S 5 coordinates of the Wilson lines are at antipodal points on opposite sides of the contour. In our setup introduced in Section 3.1, we may choosen =n 0 for s ∈ (−T /2, T /2), and n = −n 0 for s ∈ (T /2, 3T /2). This is a perfectly sensible configuration on the N = 4 SYM 5 One may wonder whether there are any mass renormalization effects induced by the self interactions of the 1/2 BPS Wilson lines at each side of the loop, which would also have to be factored out from the Wilson loop if one wants to isolate the interaction energy between the heavy quarks. However, the supersymmetric nature of N = 4 SYM sets these corrections to zero, which may be verified perturbatively [92], and hence the bare mass subtraction is equivalent to subtracting twice the physical mass of the heavy quark from the energy of the quark-string configuration. side, and it actually preserves a maximal number of supersymmetry charges [86,92,95]. Furthermore, it immediately satisfies W [C] = 1, with C the contour introduced in Eq. (3.4), and it respects the "zig-zag" symmetry, in the sense that if we extend the two timelike segments by an arbitrary extent, the contributions from each side of the contour cancel each other.
Perhaps the most intuitive argument for the time evolution of a heavy quark-antiquark pair to be represented by two Wilson lines with antipodal positions on S 5 comes from inspecting their purported equations of motion. 6 Namely, we can define the notion of a heavy antiquark as the object that transforms in the representation conjugate to that of a heavy quark, and therefore follows a evolution equation conjugate to that of the heavy quark. This means that if we take the evolution equation for a heavy quark Q (with its mass already subtracted) to be then the evolution of a heavy antiquark follows where the arrows on top of ∂ 0 indicate the directions for it to act. If one then constructs the supersymmetric Wilson loop (3.13) that describes the joint evolution of this heavy quark-antiquark pair, the sign flip in front of A 0 is accounted for by flipping the sign ofẋ µ along the same path, and the sign flip in front of ϕ is accounted by inverting the direction n. Therefore, based on these considerations, we have a good candidate for the analogous object to the QCD chromoelectric correlator (3.12) to calculate in the N = 4 SYM theory, namely, the correlation function obtained from a 1/2 BPS Wilson loop that has two timelike lines at antipodal positions on S 5 , which has desirable properties distinct from the loop with constantn. That being said, we will nonetheless evaluate the resulting correlation functions from both of them in the following sections, for two reasons: 1. Even if the limit L → 0 of the expectation value of the Wilson loop with constantn does not satisfy our expectations, it might still prove instructive to study non-Abelian electric field insertions in a Wilson loop that describes the interaction energy between a pair of heavy quarks; 2. Ultimately, the motivation behind this calculation is the phenomenology we want to extract for heavy particle pairs in a thermal medium. Since we are using this holographic setup as a model of QCD, we should evaluate all objects that have a reasonable chance to resemble our correlation function of interest, and judge them by their phenomenological implications.
Before proceeding, we also want to comment on a more recent prescription [98,99] to evaluate the standard Wilson loop (3.1), which arguably should take precedence in our analysis over the locally 1/2 BPS Wilson loop because it is constructed from exactly the same fields as in the standard Yang-Mills Wilson loop. This prescription states that the strong coupling limit of Eq. (3.1) is given by extremizing the Nambu-Goto action with Neumann boundary conditions on the S 5 . In practice, this means that most results for the standard Wilson loop at strong coupling are the same as those for the Wilson loop (3.13) with constantn, and they only start to differ at the next order in 1/ √ λ [100]. This is true because, incidentally, a constantn configuration is consistent with Neumann boundary conditions on S 5 .
However, the situation for the limit of our interest, namely, the construction of a standard Wilson loop with two overlapping anti-parallel timelike Wilson lines such that W [C] = 1, using the Neumann boundary condition prescription, might be more subtle. To see this, let us re-examine the arguments that gave rise to this prescription, and to the conclusion that the strong-coupling results of Eqs. (3.1) and (3.13) are the same for a general Wilson loop. To facilitate the references to previous works, we will work in Euclidean signature for the remainder of this subsection, unless otherwise noted or if it is explicit from the discussion (e.g., if we refer to lightlike or timelike).
The first reference to Neumann boundary conditions was given by Drukker, Gross and Ooguri in Ref. [92], where they proposed a boundary condition for an even more general Wilson loop: where y i = y i (s) ∈ R 6 is now a general vector. Concretely, the boundary condition they prescribed was where X µ denotes the usual Mink 4 coordinates and Y i is a 6-dimensional vector with the magnitude given by the AdS 5 radial coordinate z, and direction specified byn ∈ S 5 , h ab = ∂ a X µ ∂ b X ν g µν is the induced metric on the worldsheet, σ 1 is the coordinate parametrizing the boundary contour, and σ 2 is a coordinate that runs into the worldsheet. An immediate consequence of these boundary conditions, which is discussed explicitly in Ref. [92], is that the area-minimizing extremal surface reaches the boundary z = 0 of AdS 5 if and only if the loop variables obey the constraintẋ 2 =ẏ 2 . That is to say, only in this case σ 2 = 0 corresponds to z = 0. Once this constraint is incorporated, the inhomogeneous Neumann conditions in Eq. (3.22) become Dirichlet conditions on the S 5 that selectn i (σ 1 , z = 0) =ẏ i /|ẏ|.
The boundary condition first proposed in Ref. [98] for the pure gauge Wilson loop (3.1), i.e., homogeneous Neumann boundary conditions for the S 5 coordinates, is exactly of the same form as in Eq. (3.22) for the S 5 variables with the right hand side vanishing, but prescribes z = 0 as a Dirichlet condition. This is more explicitly written in Ref. [99], where the boundary conditions are written as z = 0 and n a h ab ∂ b U i = 0, where, in their notation, U i ∈ S 5 is the unit vector that we calln, and n a is a unit vector normal to the worldsheet boundary. 7 The motivation behind this prescription of homogeneous boundary conditions for the S 5 variables in Ref. [98] was not disconnected from the preceding discussion of the inhomogeneous Neumann boundary condition shown in Eq. (3.22) introduced in Ref. [92]. Indeed, part of the reasoning that led to this prescription in Ref. [98] was that when we go to the real time description of Eq. (3.20) (i.e., to Minkowski signature), provided the conditioṅ y 2 =ẋ 2 is met, the coupling to the scalars disappears if we choose x µ to be a lightlike patḣ x 2 = 0, because in this way we forceẏ i = 0, and as such, the original Wilson loop (3.1) is recovered. 8 The boundary conditions induced on the S 5 variables by this boundary contour are exactly homogeneous Neumann conditions, thus substantiating the proposal in Ref. [98]. It then also follows that, if the path of a supersymmetric Wilson loop at constant n can be approximated by a lightlike path in such a way that the extremal worldsheets of both configurations are approximately the same, converging onto each other when the lightlike approximation of the original path becomes better and better (e.g. by using many small lightlike segments with neighboring segments perpendicular to approximate a straight line), then the expectation values of both Wilson loops in Eqs. (3.1) and (3.13) calculated with their respective prescriptions will agree, since a constantn fulfills the homogeneous Neumann condition.
However, for our setup where the two anti-parallel Wilson lines are coincident in space, ifn is held constant, it is not clear whether a lightlike deformation of the (timelike) boundary contour produces a controllable approximation to the worldsheet obtained from the original undeformed contour. Indeed, given that the radial AdS 5 extent of the worldsheet we will find in Section 4 is proportional to the distance L between the two Wilson lines, which we want to take to zero L → 0, any deformation (lightlike or otherwise) will qualitatively affect the resulting extremal surface, and thus there is no guarantee that the generalized area of the worldsheet with the deformed contour will have the same value as the original undeformed one. Furthermore, for the proper definition of our observable, the first limit we should take is L → 0, because it is only in this limit where we are free to modify the temporal extent of the contour along the time direction without changing the result. This motivates us to look more closely at how the homogeneous Neumann condition can be obtained as a dual description of the standard Wilson loop (3.1) and the role ofn.
Another way of arguing for the Neumann prescription for the pure gauge Wilson loop (3.1) has been given in Ref. [99]. Namely, by considering the value ofn on the boundary contour as another dynamical variable to be integrated over in the path integral, one can write

23)
7 Choosing coordinates such that na ↔ ∂/∂σ2, one can show that nah ab ∝ h1aϵ ab , and therefore both ways of writing down the Neumann boundary condition are equivalent. 8 See Section 4 of Ref. [98].
where both ⟨W BPS [C;n(s)]⟩ and e iS NG [Σ(C;n)] are obtained from a path integral in N = 4 SYM and 10-dimensional Supergravity, respectively. On the right hand side of this equality, treatingn(s) as a dynamical variable gives the Neumann condition as an equation of motion.
On the left hand side, one can argue as in Ref. [99] that this path integral gives the pure gauge Wilson loop (3.1). This is achieved by expanding the Wilson loop in a power series inn(s), noting that the first term is exactly the Wilson loop (3.1), and the rest of the terms either vanish by symmetry or are irrelevant operators. In fact, in both ways of arguing for the Neumann prescription, no reference to the constantn solution is made in its formulation. The connection to the constantn solution only appears by noting that, since a constantn satisfies the Neumann condition, it follows that the dual descriptions of Eqs. (3.1) and (3.13) are governed by the same saddle points, provided that this saddle point is the dominant one, which is usually the case.
However, direct inspection of the left hand side of Eq. (3.23) reveals that there is another saddle point for our contour of interest C, given in Eq. (3.4). Namely, the equations of motion that extremize the left hand side of Eq. (3.23) have a solution given byn(s) =n 0 for s ∈ (−T /2, T /2) andn(s) = −n 0 for s ∈ (T /2, 3T /2), wheren 0 is a fixed direction on S 5 . This is easy to see: any first-order variation of the Wilson loop will give zero because varying it with respect to any field in it will result in an operator insertion along the loop, proportional to an SU(N c ) generator matrix T a . For the antipodaln(s) configuration that we claim as a solution, in the equation of motion for then(s) variable, the Wilson lines cancel and all that remains is proportional to Tr color [T a ] = 0. The rest of the saddle point configuration is determined by the variations of the N = 4 SYM action, which provides its standard equations of motion (i.e., the same equations as in the absence of a Wilson loop).
Furthermore, when we selectn with antipodal positions on the S 5 , because the Wilson lines are at separate, antipodal coordinates, it seems plausible that the contour C as given in Eq. (3.4) does admit a lightlike approximation that modifies the extremal surface in such a way that it converges to the unmodified solution as the approximation is made finer. This is so because, as opposed to the case of Section 4, the size of the deformation here can indeed be taken to be much smaller than the radial extent of the unperturbed worldsheet we will find in Section 5, which is (πT ) −1 . Proceeding in this way, we will find a saddle point that should 9 respect the homogeneous Neumann boundary condition (it is guaranteed to respect it if the proposal of Ref. [99], expressed through Eq. (3.23), holds), and yields a result ⟨W [C]⟩ = 1. Therefore, because the pure gauge Wilson loop (3.1) satisfies the unitarity bound ⟨W [C]⟩ ≤ 1, we would have necessarily found an extremal solution of minimal energy, and consequently, this saddle point would be the one that provides a dual description of the Wilson loop on the contour shown in Eq. (3.4). We stop short of claiming a proof of this result because an explicit verification should also provide the extremal worldsheet that is dual to the Wilson loop on the path C given by Eq. (3.4) and explicitly verify all of the boundary conditions discussed above. However, we do conjecture that the pure gauge adjoint Wilson line has the antipodaln configuration as its gravitational dual through the gauge/gravity duality.
To close this section, we note that when there is a nonvanishing spatial separation between the Wilson lines that comprise a Wilson loop, as in the case of the heavy quark interaction potential, one can indeed use the same solution as at constantn to describe the extremal worldsheet that gives the expectation value for the pure gauge Wilson loop (3.1), in accordance to the conjecture laid out in Ref. [98]. One may then worry about the unitarity bound ⟨W [C L ]⟩ ≤ 1 being violated. However, the situation here is rather similar to that of QCD: only after regularizing and renormalizing does the statement log⟨W [C L ]⟩ ∝ 1/L make sense. For concreteness, we consider QCD on a 4-dimensional Euclidean lattice, characterized by a lattice spacing a. If one considers simply the expectation value ⟨W [C L ]⟩ in the limit L → 0 at fixed a, the value will only converge to 1 if lattice artifacts are taken into account, which happen when L ≲ a. This means that to access ⟨W [C L=0 ]⟩ on the lattice, one has to take a → 0 first before taking L → 0. At finite L the expectation value ⟨W [C L ]⟩ contains both an L-dependent term and an L-independent diverging term, and they cancel at L = 0. However, the ratio ⟨W [C L 1 ]⟩/⟨W [C L 2 ]⟩ will be finite as a → 0 with L 1 , L 2 > 0 and will exhibit features of a Coulomb potential at distances small compared to the non-perturbative scale Λ −1 QCD . While it is true that ⟨W [C]⟩ = 1 for the contour shown in Eq. (3.4), its ratio with ⟨W [C L ]⟩ at L > 0 in the continuum limit is formally infinite. To give meaning to the Wilson loop at a finite spatial separation L and extract energy differences from ⟨W [C L ]⟩, renormalization and regularization is required. Taking ratios achieves this immediately. The need for extra regularization is also clear from calculations in the gravitational description where the boundary contour is modified to be lightlike, see, e.g., Eq. (3.7) in Ref. [98].

Variations of Wilson loops in AdS/CFT
Having introduced the necessary concepts to formulate the calculation of the expectation value of a Wilson loop and its path variations (i.e., functional derivatives) that provide the non-Abelian electric field correlator we want to calculate in Yang-Mills theory, we now proceed to describe how the calculation of these path variations takes place in the gravitational description of a Wilson loop in N = 4 SYM. Concretely, we will first lay out the steps one needs to follow to extract correlation functions from the calculation of an extremal worldsheet in a gravitational description and its derivatives. We will defer their evaluation to the concrete setups we discuss in sections 4 and 5.
Secondly, since we are interested in the real-time correlation functions of operators in a thermal ensemble, we will discuss the setup of our calculation on the Schwinger-Keldysh contour and how it is realized holographically in the dual gravitational description in the last part of this section. Specifically, we will discuss the iϵ prescription appropriate for time-ordered quantities, which is exactly the nature of the correlation function we want to calculate. We will also discuss the qualitative differences between the observable we will calculate and the correlation function that defines the heavy quark diffusion coefficient.

Generating non-Abelian electric field insertions
Our goal is to insert field strength operators along the boundary contour C. As we described in Section 3.1, this is achieved by taking functional derivatives with respect to deformations of the contour C, parametrized by h i (t). The corresponding operation in the gravitational description is to take functional derivatives with respect to the boundary conditions of the extremal surface. Operationally, we can achieve this by: This can be done systematically, starting from the lowest order up to the desired number of powers in the perturbation. The result may be organized as where the kernel can be obtained from the four steps listed above.
With this definition, S (p) NG [C; h] is the generating functional for (non-Abelian electric) field strength insertions along the contour C up to order p, which allows one to evaluate correlation functions with up to p insertions of operators. This object fully characterizes the n-point functions of non-Abelian electric field strength insertions along the contour C at leading order in the large 't Hooft coupling limit, as discussed in the paragraph before Eq. (3.16). To marginally ease the notation, let us introduce where C h denotes the contour C perturbed by h. We note, however, that the Wilson loop in the duality (3.14) is not the pure gauge loop (3.1), because it also involves scalar fields (3.13), and so its path variations also give rise to a contribution from the scalars that modifies the chromoelectric field operators. Namely, assuming that the fluctuations h i are only in the spatial directions, the non-Abelian field strength insertions are modified tȯ wheren is the direction on the S 5 that appears in (3.13), and [D i ϕ] a = ∂ i ϕ a + gf abc A b i ϕ c is the gauge covariant derivative of the scalar field. Therefore, after defining what we can calculate using the duality (3.14) is In this expression,n takes the sign that it has on the part of the contour where the derivative is taken. 10 If our conjecture at the end of Section 3.2.1 holds true, then we also have ∆ ij (t) = ∆ EE ij (t), and the distinction becomes idle. In this case, according to the prescription presented in Ref. [98], the contribution from the scalars disappears.
In these expressions, we have included an extra normalization factor given by the unperturbed Wilson loop. As explained before, we should have ⟨T W [C]⟩ = 1 for a Wilson loop in pure Yang-Mills theory (or QCD), and also ⟨T W BPS [C]⟩ = 1 for the configuration we discuss in Section 5, but, as we will remind the reader later, this is not the case for the loop that describes the heavy quark interaction potential, which we discuss in Section 4. In this situation it is appropriate to normalize the correlator by the expectation value of the "background" Wilson loop.
When we do have ⟨T W [C]⟩ = 1, we can summarize the above result as where we have used the standard normalization T F = 1/2 for the fundamental representation of SU(N c ). Therefore, the kernel ∆ ij is exactly the object we are interested in, for each worldsheet configuration of interest. Before proceeding to calculate them, we will take a small digression to discuss aspects of the worldsheet near the turning points at t = ±T /2. In particular, we will discuss how the time-ordering iϵ prescription emerges in this setup, and we will also comment on the interplay between the chosen form for the boundary conditions h i and how the fluctuations behave near t = ±T /2.

The Schwinger-Keldysh contour in AdS/CFT
The fact that we want to calculate a thermal expectation value requires us to introduce the Schwinger-Keldysh (SK) contour [101] in order to represent the observable of interest through path integrals. The holographic realization of the Schwinger-Keldysh contour in AdS/CFT dates back to the early work of Herzog, Son and Starinets [102,103], which was more recently expanded and refined by Skenderis and van Rees [93,94,104]. In a nutshell, each segment of the Schwinger-Keldysh contour is the boundary of an asymptotically AdS 5 bulk geometry, and these bulk spacetimes are glued together according to appropriate matching conditions, discussed at length in Ref. [93]. This allows one to formulate the Schwinger-Keldysh contour (and, consequently, the bulk geometries) in the complex time plane holographically. As in any quantum-mechanical theory, the thermal nature of the average is dictated by the fact that modes with energy ω are coupled to themselves across a Euclidean time direction of extent β, which gives rise to the characteristic thermal statistics through factors of e −βω , and this is realized in the holographic setup by matching the bulk manifolds accordingly. Importantly, the contour in the complex time plane also defines the iϵ prescriptions necessary to define the correlation functions in the limit where we take T → ∞. This limit is convenient to do calculations because it restores local translational invariance in the time direction, which is also necessary if one wants to extract the Fourier components of a correlation function at an arbitrary frequency ω, and crucially, to have a continuous limit of it as ω → 0. Because of this, and as we will state explicitly later on, we will indeed take the limit T → ∞ for both of the configurations that we have introduced and will study in Sections 4 and 5.
We will derive the appropriate iϵ prescription for our setup. Before proceeding, it is also important to note that not all observables are equally sensitive to the thermal nature of the contour. In particular, for our correlation function of interest, which is defined through a Wilson loop that "backtracks" over the path that it traverses to cover the distance between the two electric field insertions, the extremal surface that defines the expectation value of this Wilson loop will lie only inside the bulk manifold that has the time-ordered part of the Schwinger-Keldysh contour as its boundary. This is so because, crucially, all operators are time-ordered in our setup. The main consequence of this is that the fluctuations that propagate on top of this extremal worldsheet will lose sensitivity to the Euclidean (thermal) part of the Schwinger-Keldysh contour. This is one of the most important qualitative differences between our setup and the heavy quark diffusion calculation in Ref. [65], where such a distinction is necessary. We will elaborate further on this at the end of this section, after establishing the prescriptions to select the mode solutions in our setup, in accordance with the boundary conditions that the fluctuations must satisfy.
The iϵ prescription for time-ordered correlation functions The following discussion applies to both configurations to be discussed in Sections 4 and 5. We will first make use of the deformability of the Schwinger-Keldysh contour to derive the iϵ prescription that is appropriate for the time-ordered correlation function we want to calculate. Secondly, and as a consistency check, we will verify that this prescription is the same as that one would get by considering the contributions from the turnaround points (t = ±T /2) of the Wilson loop where the path C becomes spacelike.
To carry out this derivation, it is first necessary to make some remarks about the nature of the Schwinger-Keldysh contour. In principle, any contour starting at t = t i and ending at t = t i − iβ in the complex time plane is adequate in order to evaluate thermal expectation values. The only requisites are that: 1. The operators of interest are placed at some point along this contour, and 2. The contour itself never goes upward in the complex time plane, i.e., its tangent vector always has a non-positive imaginary part. This is necessary for the path integral to be evaluated by a saddle point approximation.
The standard Schwinger-Keldysh contour achieves this by first going from t = t i along the real axis to some final time t f (which defines its time-ordered segment), then turning around and going from t f to t i along the real time axis (infinitesimally displaced by −iϵ, which defines the anti-time-ordered segment), and finally going from t i to t i − iβ to close the contour. In fact, the second point of our requisites suggests that we can tilt the timeordered and the anti-time-ordered parts of the Schwinger-Keldysh contour slightly: Instead of drawing the time-ordered contour exactly along the real axis, parametrized by is the parameter running along the contour, the natural way to have a well-defined saddle point is to take t c = t i + t ′ (1 − iϵ) for the time-ordered branch, with ϵ a small (infinitesimal) positive number. The anti-time-ordered branch, going back to t i must then be parametrized by The contour is then closed by a straight line along the imaginary time axis at Re(t c ) = t i downward until reaching t = t i − iβ. As we will see momentarily, these iϵ deformations provide exactly the standard time-ordered and anti-time ordered prescriptions to evaluate correlation functions.
With this setup, we can consider the Nambu-Goto action that describes the worldsheet in the spacetime that has the time-ordered segment of the Schwinger-Keldysh contour as its boundary. In terms of the parameter t ′ , which we write as t in what follows, the metric reads as The more general case we will consider for our purposes is that of a worldsheet parametrized by X µ = (t(s, z), x(s, z), y(s, z), 0, z,n(s, z))) , (3.32) where s is a worldsheet coordinate that parametrizes the Wilson loop, which we may define at z = 0 to be the arc length on the loop. The spacetime coordinates t, x, z,n describe the background solution, while y describes the fluctuations we seek to solve for. In our setup, n is completely determined by a large circle angle ϕ(s, z) due to symmetry considerations.
In what follows, we will denote the derivatives of a coordinate a by da dz = a ′ and da ds =ȧ. Since we will be considering essentially infinitesimal perturbations y, whether the worldsheet is spacelike or timelike is wholly determined by the background solution. Expanding up to quadratic order on y, one obtains the following action: Moreover, Eq. (3.35) provides a prescription that affects the equations of motion of y. To see this in a concrete setup (which will be relevant in Section 5), we consider the case of a radially infalling background worldsheet at constant x andn, with t = s. The action for the fluctuations simplifies to 36) which implies that at the level of the equations of motion, the frequency ω of the mode solutions will always appear as ω 2 (1 + iϵ). This in turn defines the pole prescription to evaluate the propagator. As we will see later in Section 5.3, this also determines which mode solution should be used when calculating correlation functions. One may also wonder whether the behavior of the worldsheet around the turnaround times t = ±T /2 affects this conclusion. Specifically, one can wonder whether one can get extra imaginary terms in the equations of motion by having a transition where the induced metric on the worldsheet goes from having Minkowski signature (i.e., timelike) to having Euclidean signature (i.e., spacelike).
To have control over the behavior of the worldsheet at the turnaround times t = ±T /2, we need to regulate the backtracking of the loop in a way that its tangent vector is continuous throughout, such that if we look closely enough, the extremal surface will still be smooth. Our choice of regulator for the present purpose is to introduce a small spatial separation L between the two lines, which is compatible with the discussion in the previous sections. Conversely, the only way to smoothly turn from a timelike tangent vectorẋ µ going in the future direction to one going in the past direction is by having a segment where it is spacelike. As such, our choice for a regulator is actually generic.
This motivates studying the behavior of a worldsheet close to a spacelike boundary segment. To gain intuition, let us first discuss a few examples. A family of solutions that is easily obtained at These solutions are spacelike surfaces that satisfy the Euler-Lagrange equations obtained from the Nambu-Goto action that are bounded by the hyperbola t 2 − x 2 = ρ 2 0 at z = 0. Another family of solutions to the Euler-Lagrange equations is given (implicitly) by the integral where z c is a parameter defining different solutions, all of which are bounded by the line t = 0 at z = 0, valid in a small neighborhood of a spatial Wilson line segment with varying x and all else held constant (to fully determine a unique solution it is necessary to specify how the surface is closed, or equivalently, how the Wilson loop path closes itself, far away from the region we just studied). All of these have the crucial property that they are spacelike surfaces, which motivates investigating whether our previous conclusion is affected when we deform the contour slightly by introducing a spatial separation.
Note that the iϵ in Eq. (3.34) also provides a prescription to evaluate the action in the case of a spacelike worldsheet. Specifically, it determines that a spacelike worldsheet has a Nambu-Goto action determined by the substitution which, satisfactorily, is exactly what we would get by demanding that whenever the worldsheet is spacelike, the Nambu-Goto action should be the same as if we had started in Euclidean signature from the beginning. Now we may ask what happens if we include perturbations on top of a background worldsheet that features a transition from spacelike to timelike and vice-versa. Given that these perturbations are introduced on top of a solution that extremizes the action, the action for the fluctuations in a spacelike region should be real and positive definite. We will verify this explicitly in what follows, as it will be crucial to our results that the iϵ prescription would not be modified by contributions from a spacelike region.
When the background worldsheet is spacelike, the argument of the square root in Eq. (3.34) becomes negative, and we must therefore use Eq. (3.37) to get which for the fluctuations read (3.39) As written, this is a general expression. However, the expression is explicit enough for us to make generic statements about how the fluctuations y(s, z) behave on a background specified by X µ = (t(s, z), x(s, z), 0, 0, z,n(s, z)). The key observation is that the quadratic form in the numerator of the integrand in Eq. (3.39) can be written as 40) and noting that the 2 × 2 matrix in the middle of this expression is equal, component , 0, 0, z,n(s, z)) describes the background solution, with the first component of the matrix (for the indices α, β) corresponding to a derivative with respect to s, and the second with respect to z. Therefore, if the background worldsheet is spacelike, it follows that both eigenvalues of the induced metric ∂ α X µ ∂ β X ν g µν are positive, and hence, that it is a positive definite matrix. Consequently, the action for the fluctuations (3.39) is positive definite whenever the background worldsheet is spacelike. Then, extending the boundary contour as T → ∞, the contributions of these regions will be of the form (for definiteness, consider for some positive definite quadratic form Σ (note that the minus sign in the definition (3.33) means that the Gaussian integral over the fluctuations y is convergent). The positive definiteness of Σ is guaranteed by the fact that it is determined by the induced metric of a spacelike surface, as discussed below Eq. (3.39). The same is true for the region at τ = −∞.
The net effect of this on the action, after decomposing y in terms of the mode functions at intermediate times, is to add an iϵ to the ω 2 coming from the temporal derivatives in the action, with ϵ > 0. 11 This iϵ modifies the mode equations by effectively shifting ω 2 → ω 2 (1 + iϵ), in the same way as the modification induced by the Schwinger-Keldysh contour tilts on the complex plane. As such, the prescription is unambiguously determined. As a side note, we remark that the above discussion did not address how the solutions on the different sides of the turnaround region are coupled to each other. Continuity of the fluctuations and of their (appropriately normalized) derivatives is a first requirement, but, as hinted from our previous discussions, closer inspection from the field theory perspective reveals that the solutions must actually be more constrained than that. To see this, if instead of considering perturbations as given by Eq. (3.8), we consider we then obtain W [C f ] = 1 for all f , since the loop consists of a single line that was traveled back and forth on top of each other. Hence, taking derivatives with respect to this kind of contour deformations gives zero, and as such, the corresponding response kernel for Wilson loop variations evaluated with AdS/CFT techniques must also be identically zero, whenever the loop satisfies W [C f ] = 1 to begin with. We stress that this is the case for the loop witĥ n at antipodal positions on the S 5 , but it will not necessarily be the case whenn is constant (even though, as we will see later, deformations as in Eq. Differences with the heavy quark diffusion coefficient Finally, let us comment on how this calculation differs from the one for the heavy quark diffusion coefficient [65]. To do this, it is most helpful to use the construction put forth by Skenderis and van Rees [93,94,104], where the Schwinger-Keldysh contour has a concrete holographic realization. This is done by constructing a bulk manifold made up of several submanifolds satisfying appropriate matching conditions, where the boundary of each of these submanifolds is identified with the lower-dimensional spacetime that corresponds to a given segment of the Schwinger-Keldysh contour in the boundary CFT. This can be done in the same way for the correlator we are presently considering and the one that determines the heavy quark diffusion coefficient. However, the manifold on which the fluctuations that are of interest to us propagate is not the full manifold associated to the Schwinger-Keldysh contour of the full CFT. Rather, the fluctuations propagate on a lower-dimensional manifold given by the background solution for the worldsheet configuration. Whether this worldsheet configuration spans every region of the Schwinger-Keldysh contour is solely determined by the shape of the boundary Wilson loop. In the case of the heavy quark diffusion coefficient setup, the corresponding Wilson loop consists of a single Wilson line that winds around the Schwinger-Keldysh contour once. This means that the background worldsheet configuration that defines the manifold on which fluctuations will propagate spans the whole bulk space at a single fixed position coordinate on the boundary x, and is parametrized by a radial AdS coordinate and a temporal coordinate that goes over both Minkowski and Euclidean regions of the bulk geometry. The topology of the boundary manifold is that of a circle, and the end points of the real-time segments have to be matched with those of the imaginary-time segments, in consistency with the Schwinger-Keldysh contour (see Fig. 1). Then, if one seeks for mode solutions for the fluctuations on top of this 2-dimensional geometry, the matching conditions discussed in Refs. [93,94,104] imply that Fourier modes e −iωt F ω (z) on the realtime segments (where F ω (z) is the radial AdS profile of a solution with a frequency ω on the boundary) have to be matched with solutions of the form e ±βω F ω (z) on the Euclidean The difference, qualitatively winding around the Schwinger-Keldysh contour segments. Therefore, factors of e βω naturally appear in the response functions. This is what gives rise to the KMS relations between the different types of correlation functions that can be calculated by introducing fluctuations and evaluating the response functions on different segments of the SK contour. 12 On the other hand, the Wilson loop that defines the correlation function we are presently interested in does not wind around the SK contour. That means that the "thermal" contributions e βω that come from matching the fluctuations around the contour will not be present in this case, and therefore all of the temperature dependence that will appear is going to be due to temperature effects on the bulk geometry of the Minkowski part of the manifold that holographically realizes the path integral associated to the Schwinger-Keldysh contour. Hence, the way in which both observables are defined has manifestly distinct effects in the values that the correlation functions take. After discussing the calculations in detail throughout the following sections, we will see how these differences manifest themselves in the final result.

The Wilson loop with constant S 5 coordinate
As we discussed in Section 3.2.1, the standard choice to do calculations of Wilson loops using the AdS/CFT correspondence in strongly coupled N = 4 SYM is to set a constant value forn throughout the Wilson loop. This is indeed the setup used in the celebrated paper by Maldacena [86] to calculate the heavy quark interaction potential at strong coupling. Since our interest is to describe the dynamics of a heavy quark-antiquark pair close together, we find this is a natural starting point that warrants exploration, regardless of our previous observation that this choice forn does not preserve all properties we expect from the Yang-Mills Wilson loop. We will consider the situation where then coordinates are at antipodal points on the S 5 for the heavy quark and the heavy antiquark respectively in Section 5.
The calculation consists of three steps. First, in Section 4.1 we will discuss the "background" worldsheet solution that hangs from the unperturbed Wilson loop (i.e., without the deformations that give rise to the field strength insertions), thus establishing the geometry on which fluctuations can propagate. Secondly, in Section 4.2 we will discuss the action for the perturbations on top of this background solution, and how to extract the correlation function of interest for the specific geometry we describe in the first step. Finally, in Section 4.3 we will calculate the correlation function as prescribed by the previous steps. We will provide more details for fluctuations that are transverse to the worldsheet, and discuss longitudinal fluctuations (to be defined in what follows) in a more succinct way. We will also check our results by a numerical calculation of the background extremal surface and the correlation function in Euclidean signature.

Background
The heavy quark interaction potential can be extracted from a rectangular Wilson loop of temporal extent T and spatial separation L, with T ≫ L. Its calculation using supersymmetric Wilson loops has been discussed many times in the literature. The original papers discussed this at zero temperature [86,87]. More general setups were later discussed including finite temperature effects and a relative velocity between the heavy-quark pair and the medium [62,[107][108][109][110][111]. The same loop has also been considered in AdS/QCD to calculate the characteristic correlation lengths of field strength correlators [112] in the limit L ≫ T . In what follows, we review the main features of the extremal surface that appears in the AdS/CFT calculation of the static heavy-quark potential in N = 4 SYM. Our goal is to study it in the limit L → 0, where the two parallel Wilson lines that construct the timelike segments of the loop get pulled close together. Because the solution for this Wilson loop has been well-studied in the literature, we will only briefly review the results and highlight their most important features for our purposes.
In the presence of a black hole described through the metric (3.16), the minimal area worldsheet configuration that hangs from a rectangular contour of size T ×L on the boundary, with T ≫ L, can be parametrized by where σ ∈ [0, L] and τ ∈ [−T , T ]. For such a parametrization, the Nambu-Goto action reads Because the action does not depend on σ explicitly, there is a conserved quantity, which is the analog of the Hamiltonian H in standard classical mechanics, with H = pq − L and p = ∂L/∂q. Using this conserved quantity, one finds that the background worldsheet satisfies where we have denoted f m = f (z m ). This equation can be integrated to find an implicit solution for z(σ), which is given by where z m is the maximum value of the radial AdS coordinate z(σ) and F 1 is the Appell hypergeometric function. It is in turn given by This equation has two solutions for any given value of πT L < πT L max ≈ 0.86912, corresponding to a value of z m given by z crit m ≈ 0.84978. This is depicted in Fig. 2. As discussed in Refs. [113,114], the solutions with z m > z crit m are unstable, and beyond this value the preferred configuration is that of two disconnected, radially infalling surfaces from two parallel Wilson lines. Since we will be interested in the L → 0 limit, namely, LπT ≪ 1, the solution we have to consider is always in the branch with z m < z crit m . The energy of this configuration 13 is given by This is a Coulomb-like potential, which diverges in the short-distance limit L → 0. The constant term proportional to T comes from subtracting the area of the disconnected worldsheet that hangs only down to the horizon z = (πT ) −1 , instead of all the way to z → ∞. This term √ λT corresponds to twice the thermal correction to heavy quark mass. The above equation means that, in the absence of another extra normalization factor in the RHS of Eq. (3.12), the result of taking the limit L → 0 will be ill-defined. As such, to  With the background solution in hand, we can now consider perturbations on top of it. By evaluating the second derivative of the action with respect to the perturbations, we can extract the non-Abelian electric field correlation function we seek.

Fluctuations
We now describe the dynamics induced on the worldsheet by small deformations on the contour that bounds it. Following the discussion we presented in Section 3.3.1, one arrives at the conclusion that this is achieved by introducing small fluctuation fields on the string, which capture how the boundary perturbations propagate into the string. These fields obey second-order partial differential equations that are determined from the Nambu-Goto action. Once one solves the corresponding differential equations, one has to evaluate the Nambu-Goto action "on-shell," i.e., on the solution to the equations of motion, as a function of a general boundary condition h i (t). Then, by taking functional derivatives with respect to h i (t) of the on-shell action, one can extract the correlation function we are interested in. Because of notational clarity, we will give the results in terms of the kernel ∆ ij , which will be different for this configuration than that for the configuration in Section 5. We will denote this section's expression for this kernel by ∆ c ij , where the "c" may stand for "connected" or "Coulomb," in reference to the shape of the background worldsheet and to the nature of the interaction potential, respectively. We will denote the solution of the next section by ∆ d ij with "d" standing for "disconnected". Only after we have both of them at hand will we use Eq. (3.30) to relate our answer to the non-Abelian electric field correlator of interest.
To evaluate ∆ c ij , the first step to take is introduce perturbations along all of the AdS 5 coordinates, and consider a string parametrized by X µ (τ, σ) = (τ + y 0 (τ, σ), σ + y 1 (τ, σ), y 2 (τ, σ), y 3 (τ, σ), z(σ) + y 4 (τ, σ),n) . (4.7) We do not consider fluctuations on the S 5 coordinates because they are decoupled from the rest at the quadratic level in the Nambu-Goto action, which is all we need to evaluate our correlator. However, this parametrization has redundancies in it, because fluctuations that lie on the tangent space to the worldsheet are not physical deformations, but rather a coordinate reparametrization. This means we can choose our worldsheet coordinates such that we can set y 0 (τ, σ) = 0, as well as set to zero a certain linear combination of y 1 and y 4 . To find it, we need to project y 1 and y 4 along the directions parallel and perpendicular to the worldsheet. Let δ(τ, σ) parametrize the fluctuations orthogonal to the worldsheet, and r(τ, σ) describe reparametrizations along the worldsheet. Projecting along the parallel and orthogonal directions to the tangent vector of the worldsheet by means of the metric g µν , one finds that the parallel and perpendicular fluctuations are parametrized by As long as we are in the linear response regime, it is a straightforward exercise to show that the Nambu-Goto action only depends on δ ⊥ X µ , with no dependence on δ ∥ X µ after using the background equations of motion. 14 Therefore, we can describe the perturbed worldsheet by ,n) . (4.10) This will be accurate as long as the fluctuations can be treated perturbatively, which is indeed the case of interest because we only need to evaluate the derivative of the action about the background configuration up to quadratic order in the fluctuations, which means that the corresponding equations of motion we will have to solve are linear. The next step is to write down the action up to quadratic order and derive the equations of motion for the perturbations. After using the background equations of motion and the boundary conditions, one finds that the linear terms in the fluctuations vanish, and it is then straightforward to show Eq. (3.15) becomes where each term is given by (4.14) In these equations, ∥ and ⊥ should be distinguished from the meanings of being tangent and perpendicular to the background worldsheet. Rather, they indicate whether the perturbations on the boundary (z = 0) are in the same plane as the Wilson loop or perpendicular to it. From the action S NG,c , one can derive the equations of motion for the fluctuations. We want to emphasize that if we had kept the redundant fluctuations y 0 (τ, σ), r(τ, σ), we would have obtained an action containing them up to quadratic order. However, upon calculating their equations of motion, one finds that they are trivial (they vanish identically), and also do not enter the equations of motion for the rest of the fluctuations up to the linear response level.
Before proceeding to the calculation of the kernel ∆ c ij , we note that from here we can already give formal expressions for the on-shell action at quadratic order in the perturbations. As one can always do for a quadratic action of dynamical variables and their first derivatives, we can integrate the action density by parts to obtain the equation of motion plus a total derivative, which reduces to a boundary term. Using Eq. (4.3) and considering nonzero boundary conditions at the timelike segments of the Wilson loop only, we find (in the limit T → ∞) (4.16) which conveniently are of the same form. The relative simplification of these expressions is due to the fact that all of the coefficients of δ and δ ′ are evaluated at the boundary z = 0.
The remaining task is to find the derivatives δ ′ and y ′ that solve the equations of motion derived from Eqs. (4.13) and (4.14), in terms of general boundary conditions on the timelike segments of the Wilson loop, of the form implied by Eq. (3.8). Concretely, we seek δ ′ and y ′ whose boundary conditions at the timelike segments of the Wilson loop are given by  There is no sign flip in the boundary condition for δ relative to that of y because the parametrization (4.10) already takes it into account.
Because the corresponding equations of motion are linear, it is possible to write down the derivative of the solutions at the boundaries in terms of linear response kernels K c ∥ (τ, τ ′ ), K c ⊥ (τ, τ ′ ). Because of how we have parametrized the longitudinal fluctuations, δ will be an even function of σ around σ = L/2, and y will be odd. Then, we can write with which the on-shell actions can be written as Because of the time translational symmetry in the limit T → ∞, we also have K c . From here, it is clear that by evaluating K c ∥ , K c ⊥ we will have all the information we need to evaluate the contribution of each type of fluctuations to ∆ c ij : and for the response kernels by where is the Fourier representation of the response kernel for either type of perturbation. With this, we can simply write down 29) and the problem is reduced to finding the respective response function K c ⊥/∥ (ω) at each frequency.

Calculation of the time-ordered non-Abelian electric field correlator
After setting up all of the machinery, we now describe the calculation of the response kernels for fluctuations in the configuration where the two timelike segments of the SYM Wilson loop have the same S 5 coordinates. We proceed with a greater level of detail for transverse fluctuations, which is arguably the simpler case, in the hope that it will make the longitudinal discussion less cumbersome. We also provide an increased level of detail in the hope that future calculations of fluctuations on top of extremal worldsheets to extract correlation functions from holography may benefit from this discussion.
Furthermore, in anticipation of obtaining results that might require careful regularization, we will also carry out the calculation allowing for more flexibility in the fluctuations than using a single perturbation h i . Specifically, we will set boundary conditions on the two timelike segments of the contour C independently. The purpose of this will be to verify that the contributions to the response kernel proportional to a Dirac delta function in time are not due to (omitted) contact terms in the RHS of Eq. (3.5).

Transverse fluctuations
As we just discussed, our goal now is to solve for y ω (σ) and extract its derivatives on the boundary. Varying S (2),⊥ NG,c with respect to y and transforming to the frequency domain, the equation we have to solve is where z(σ) is determined by solving subject to z ′ = 0 ⇐⇒ z = z m and z(σ = L) = z(σ = 0) = 0. In the interval σ ∈ (0, L/2) we take the plus sign (as the worldsheet goes into AdS 5 ), and the minus sign when σ ∈ (L/2, L). The most useful form of the above expression is To avoid introducing unnecessary numerical uncertainties, the best alternative is to transform the equation for y ω (σ) into an equation for y ω (z), because then we will not need to solve for z(σ) explicitly. 15 Performing the transformation, we have At this point, it is useful to introduce a rescaling of the radial AdS coordinate: ξ = z/z m ∈ (0, 1). In terms of this variable, we have The same equation applies for both copies of the transformed intervals σ ∈ (0, L/2) and σ ∈ (L/2, L). Let us denote the corresponding solutions as a function of z by y L ω (z) and y R ω (z), respectively, where L, R stand for "left" and "right" sides of the worldsheet. All that we need to specify in order to close the system are the boundary conditions. In terms of the original coordinate σ, we have which, in terms of the z coordinate, transform to To implement these matching conditions, the best way is to do a WKB-type analysis to extract the leading/possibly singular behavior of the solution near the horizon. We can do that analytically by writing where ∂F ± ω /∂ξ is finite at the turning point ξ = 1. With this decomposition, the matching conditions translate into To solve the system in terms of one of the amplitudes, we need to specify the boundary conditions. To extract the correlation function we are interested in, the natural choice is to prescribe y L ω (ξ = 0) = −y R ω (ξ = 0). However, in order to illustrate the nature of contact divergences that will appear in this calculation, we will instead consider the boundary condition y R ω (ξ = 0) = 0. We can obtain the boundary condition that defines our correlation function (i.e., y L ω (ξ = 0) = −y R ω (ξ = 0)) by taking linear superpositions of the boundary condition y L/R ω (ξ = 0) = 0 and using that the equations we are looking at are symmetric under the exchange of the L, R labels. With this, it is appropriate to define the response functions K AB (ω) as the derivative responses y ′ ω (σ) on side A due to a unit perturbation on side B.
Then, the null boundary condition y R ω (ξ = 0) = 0 at σ = L translates into which, together with the matching conditions at ξ = 1, fully determine the solution up to an overall constant: Formally, all that remains is to evaluate the response kernels K RL,c ⊥ and K LL,c ⊥ . We remind the reader that the superscripts are there to make it explicit that they represent partial contributions to the response kernel we want to evaluate, each coming from specific boundary conditions and response locations. These are determined by Furthermore, choosing the normalization of the mode functions such that y ± ω (ξ = 0) = 1, and writing the result in terms of the regular functions F ± ω (ξ) whenever possible, we have the response kernels in each case, before subtractions and regularizations, given by . (4.49) Note that all of the above expressions are valid for arbitrary L > 0, and furthermore, all of the discussion in this section holds for arbitrary, complex ω. 16 No approximations have been made. All that remains is to solve the equation for the modes y ± ω (ξ), or equivalently F ± ω (ξ), and with that the above expressions can be calculated explicitly.
It is useful to note that the denominators in the expressions for K AB,c ⊥ can be written in terms of the Wronskian of the differential equation for y ± ω . First we observe that where (4.51) The constant C can be fixed by looking at the behavior of y ± ω when ξ → 1. The result is where we have worked under the normalization y ± ω (ξ = 0) = 1. Finally, one can integrate the equation for the Wronskian to derive that Also, note that by construction we have y + ω (ξ)y − ω (ξ) = F + ω (ξ)F − ω (ξ), and all contributions to the response kernels K RL/LL,c ⊥ can be written entirely in terms of the regular functions F ± ω . We then find the denominators in the expressions for K ⊥ are given by: This is consequential because it provides a clean expression to implement the time-ordering prescription to evaluate the correlator. Concretely, the time-ordering prescription is implemented by taking ω → ω(1 + iϵ). It is then convenient to define ϕ ω (z m ) as with which the iϵ prescription implies that we can write 17 (4.56) Using this, and the fact that the mode functions F ± ω satisfy ∂ 2 F ± ω ∂ξ 2 (ξ = 0) = ω 2 z 2 m F ± ω (ξ = 0), an equality that follows from Eq. (4.34), we get where we have used our expression for the Wronskian as given by Eqs. (4.51) and (4.52).
The Wronskian is what allowed us to write the difference of the derivatives of F + ω and F − ω purely in terms of F ± ω with no derivatives. We clearly see that K RL,c ⊥ , K LL,c ⊥ are different functions. But this is fine, since we only expect them to be equal in the limit L → 0, and up to contact terms. Indeed, the last term in the expression for K LL,c ⊥ is a divergent term that is exactly of this form. On the other hand, by construction, K RL,c ⊥ will feature no such contact term contributions, because the variations of the Wilson loop, in the language of Section 3.2, are always at different values of the parameter s. While this means that the RL setup to extract the correlator gives a cleaner signal than the LL kernel, where no subtraction for contact terms is required, we shall still calculate both as a consistency check. In what follows, since we have isolated its origin, we will omit the contact term ω 2 z 2 m √ fm lim z→0 All that remains now is to evaluate the mode functions and substitute the result into the expressions for the response kernels. The equation for the mode functions F ± ω can be 17 Strictly speaking, one also has to analyze the mode functions and determine explicitly whether F + ω F − ω gives another O(ϵ) contribution when introducing the prescription. As it turns out, this contribution adds up with the naive one, giving the same overall effect.
found by explicitly substituting y ± ω = exp(∓i(· · · ))F ± ω into the equation of motion for y, given by Eq. (4.34). To optimize the notation, we introduce h ≡ πT z m and Ω ≡ ω/(πT ). The equation for F ± ω then reads Now we can proceed to study the behavior of the solutions. The defining condition we have to impose is the regularity of ∂F ± ω /∂ξ at the turning point, which, in terms of ∂ 2 F ± ω /∂ξ 2 being finite as ξ → 1 requires ∂F ± ω (ξ = 1)/∂ξ = 0, which can be seen by analyzing the most divergent pieces of Eq. (4.59) when ξ → 1. The near-boundary behavior of the mode functions also requires which can be obtained by expanding Eq. (4.59) in a power series in ξ near ξ = 0. While the regularity condition at the midpoint is in principle enough to determine the mode function up to an overall normalization, these near-boundary conditions may also be used in a numerical solution of Eq. (4.59). Now we want to take the limit L → 0. The process to do so is written out in full detail in Appendix C. We obtain that where the dominant contribution in the limit L → 0 is determined by (4.62) We have kept outside the definition of c 3 an overall factor of z 2 m √ fm (which does depend on L) for convenience to translate to the result for ∆ c ij , which has the inverse of this factor in the front (see Eq. (4.23) for comparison). The subleading 1/L contribution is different for each case (see Appendix C), i.e., c RL 1 ̸ = c LL 1 , but this presents no issue because we anyways do not expect the results to agree beyond the leading term as a function of L. On the other hand, the leading contribution is the same for both procedures, and it does not receive contributions from contact terms, because the construction of the RL kernel explicitly prevents this.
Having done the above, one finds that when we introduce anti-symmetric perturbations h ⊥ (t) as discussed in Section 3.2, the sum of the response kernels gives a total of and therefore the kernel that determines the two-point function for transverse deformations is given by The main feature of this result is that it diverges as L −3 when L → 0. This concludes our calculation of the response functions that determine the linear response of the Nambu-Goto action to transverse deformations on the boundary countour that defines the Wilson loop, for the background configuration that describes a Coulombtype potential between the heavy quarks. To complete the result, we now move on to calculate the longitudinal one, following the same steps.

Longitudinal fluctuations
Having gone through all the machinery in the previous section, we shall give an abbreviated discussion of the calculation for the case of longitudinal fluctuations. First, we note that to obtain the leading behavior that we got in the previous section, it is sufficient to work in the T = 0 case. This is clear by looking at how T appears in the solution for z(σ) and in the action for the fluctuations δ(τ, σ): After factoring out the overall scale L from z(σ), it is manifest that the leading appearance of T is of the order (πT L) 4 . It is then clear that, if we find a 1/L 3 dependence for ∆ c 11 (ω, L) in vacuum, this will be the dominant contribution in the limit L → 0 (note that ∆ has mass dimension three).
When T = 0, the action for the fluctuations reduces to Furthermore, if we are only after finding the leading behavior of ∆ AB,c 11 , we can even drop the term with time derivatives in this action, because we will be in the regime ω 2 L 2 ≪ 1. This will only modify the result by terms that go as 1/L. After using the conservation equation for the background worldsheet in vacuum (i.e., the conserved quantity that appears due to there not being any explicit σ dependence in the action), which we can write as z 2 √ z ′2 + 1 = z 2 m , one obtains the following equation of motion for the fluctuations: As with the transverse fluctuations, we can change variables from σ to ξ = z/z m , and rewrite this equation of motion in terms of two domains, one for σ ∈ (0, L/2), where we will use δ L , and the other for σ ∈ (L/2, L), where we will use δ R . The equation of motion for both of them is Conveniently, the solutions to Eq. (4.67) can be found explicitly: Then, the matching conditions set We can then define separate response kernels on either side for perturbations on a given side. Following our discussion of transverse fluctuations, we may define, setting δ L (ξ = 0) = A L = 1 and δ R (ξ = 0) = A R = 0, where in taking the limit we have used that the mode solutions have vanishing first and second derivatives at ξ = 0 (this can be seen directly from the mode functions, as their dependence on ξ starts at order ξ 3 ). The result is easily found to be As in the case for transverse fluctuations, these two one-sided response kernels have an equal leading order contribution to the ∆ c 11 kernel, and no contact term appears at this order. The symmetrized response kernel K c ∥ is then given by 74) which means that the contribution to the two-point function coming from the linearized fluctuations of the Nambu-Goto action is With this, we have calculated the leading contribution as L → 0 of the longitudinal fluctuations. We can therefore write the complete leading contribution to the two-point function associated to fluctuations on the extremal worldsheet that gives a Coulomb interaction potential between two heavy quarks: This completes the calculation for the quadratic fluctuations in this background configuration, and it is all we need in order to compare with the result of the next section. However, because this is highly singular as L → 0, we shall perform a numerical check that our result is not exclusive to the large T limit, and that the same behavior is obtained in the L → 0 limit for a bounded rectangular Wilson loop at fixed T .

Euclidean numerical calculation for variations on a bounded rectangle
In what follows, we will verify that the above results continue to hold when the temporal extent of the loop is finite. This is not in vain, as when T > 0 the Euclidean time direction is finite in extent, and therefore it is relevant to study the expectation value of a Wilson loop with finite temporal extentT E , to assess definitively whether the temperature can play a role in the expectation value of interest.
To demonstrate this behavior, we calculate the derivative response to transverse perturbations, solving the analogous problem to that in Section 4.3.1, but in Euclidean signature. The background solution on which the perturbations propagate is specified by the action principle where the dot stands for a derivative with respect to the rescaled imaginary timeτ E , the prime stands for a derivative with respect to the rescaled spatial coordinatex, and we have We obtain the numerical solutions to Eq. (4.77) using the pseudospectral method [115]. The pseudospectral method is an elegant way to solve boundary value problems such as the one in Eq. (4.77), which approximates the continuous differential equation by a finite set of coupled equations. Specifically, we introduce an ansatz for ξ in terms of the Chebyshev polynomials: with c ij unknown coefficients for which we need to solve and T i denoting the Chebyshev polynomial or order i. We then plug the ansatz into the equations of motion obtained from the action shown in Eq. (4.77). By evaluating these equations of motion at a finite number of points called collocation points, we obtain a set of coupled equations. The number of collocation points is given by N coll . Any collocation points that lie on the boundary of the problem are constrained using the boundary conditions instead. This immediately highlights a significant advantage of this pseudospectral method, as in this way boundary conditions are automatically satisfied, which is otherwise nontrivial for other approaches to boundary value problems.
For a general choice of such collocation points, the solution obtained in this way will not converge to the solution of the differential equation as we take the number of collocation points N coll to infinity. If, however, we choose our collocation points to lie on the simultaneous zeroes of T N coll +1 (2x − 1) and T N coll +1 (2τ E − 1), the solution obtained is guaranteed to converge to the solution of Eq. with c some positive constant. In practice, this means that with this choice of collocation points, the convergence as we take N coll → ∞ is very fast, so that we can achieve impressive precision even with a relatively small number of collocation points. By varying the number of collocation points, we can also get an estimate of our truncation error.
For linear problems, the procedure described above leads to a set of (N coll +1) 2 coupled linear equations, which is exactly the number of unknown coefficients c ij we have, so in this case one can find the solution by matrix inversion. Here we should note that the matrix that needs to be inverted is often close to singular, requiring us to work with more significant figures than machine precision provides.
A second complication is that the problem defined by Eq. (4.77) is not linear. Because of this, to find a solution we linearize the equations around a trial solution, and then use the Newton-Raphson method to iteratively update the trial solution until our iteration converges. This introduces the usual difficulties associated with Newton-Raphson, namely that for certain choices of initial trial solution the iteration might not converge, but if for the first couple of steps in the iteration one uses very small step size in the update, it is generally not hard to reach the correct solution from a reasonably chosen initial trial solution.
A sample solution for a = 0.1 and b = 10 −3 can be found in Fig. 3. We can see that for τ E away from the boundary conditions at τ E = 0 and τ E = T E , the solution agrees with the 1D problem from Eq. (4.4).
Once we obtain the background solution, we can introduce perturbations on the boundary and solve for the response functions. On each of the background surfaces, we will solve for the linear response on one timelike side of the rectangle to perturbations in the boundary conditions on the other side. We use the same decomposition as for the background solution and write the perturbation as where d ij are the unknown coefficients we need to solve for. Given that the numerical method to solve for the background worldsheet already defines a preferred basis on which to formulate this problem, we will calculate the response functions by mapping a perturbation onto a given basis element T n (2τ E − 1), where T n is a Chebyshev polynomial determining the boundary condition for the transverse fluctuations along the time axis on one side of the rectangular contour, to another basis element T m (2τ E − 1) on the other side of the contour. That is to say, given a boundary condition at x = 0, specified by y(0,τ E ) = T n (2τ E − 1), we will want to determine the response at the other side of the contour y ′ (x = 1,τ E ), decomposed in terms of Chebyshev polynomials.
To put this on a concrete mathematical footing, we shall repeat some of the discussions in Section 4.2, keeping in mind that we now work in Euclidean signature with a finite Euclidean time extent. As before, the response kernel of interest is obtained by solving the equations of motion derived from the action for the fluctuations, which can be written as .

(4.79)
Evaluating this action on-shell, with the only non-vanishing boundary conditions being in the temporal segments of the Wilson loop, we get where we have dropped ξ andξ as they vanish at the boundaries defined by the temporal segments of the Wilson loop, where by definition ξ = 0 and thereforeξ = 0. A posteriori, knowing the solution to the background worldsheet, one can verify that the limits ℓ(τ E ; a, b) ≡ limx →0,1 ξ 2 1 + ξ ′2 b 2 are finite and equal. To avoid issues with contact terms, we will extract the quadratic kernel from variations on opposite sides of the contour. That is to say, we will calculate the quadratic kernel ∆ c ⊥,E (τ E1 , τ E2 ; L) that appears as a π Γ(7/4) 3 Γ(5/4) 2 ℓ(τ E ) so that the limit is rescaled to unity.
contribution to the on-shell action, given by where K RL,c ⊥,E (τ E2 , τ E1 ) is defined as the following response kernel: All that remains is to evaluate the response kernel K RL,c ⊥,E and the limit ℓ(τ E ) in the background solution. Before proceeding, we first discuss what the expected result is. From our analysis of the effective 1D problem (in the limit T → ∞) in Minkowski signature, we see from Eq. (4.76) that the limit L ∝ b → 0 should give us Noting that δ(τ 1 − τ 2 ) = T −1 E δ(τ E1 −τ E2 ), and that in the strict limit T → ∞ we have (4.84) Given these expectations, our numerical evaluation of K RL,c ⊥,E and ℓ(τ E ) needs only to demonstrate Eq. (4.84).
We show results for b 2 /ℓ(τ E ), normalized by its limiting value in Fig. 4. We observe two general trends: 1. As b → 0, the solution indeed approaches the limiting value, converging first in the middle regionτ E ∼ 1/2 and later near the corners. For sufficiently small b, the convergence is more strongly dependent on the ratio b/a than on b alone.
2. The oscillations in the solution, which are artifacts of a truncated basis, get suppressed as we increase the number of collocation points N coll , and the convergence of b 2 /ℓ(τ E ) as b → 0 is observed even more clearly at large N coll . By increasing the number of collocation points we would reduce the truncation effects, and the solution would approach the exact profile at each value of b, with which the limit b → 0 could be examined even more precisely. However, we will refrain to go further, as we deem the plots in Fig. 4 as sufficient evidence for the value of the limit we wished to verify.
Next we test the convergence of K RL,c ⊥,E by studying the derivative response y ′ atx = 1 to a boundary condition specified by a Chebyshev polynomial T n (2τ E − 1) atx = 0. Concretely, we expand the derivative response in terms of Chebyshev polynomials and find and we plot these coefficients for different values of a and b. N pols is the number of Chebyshev polynomials we use in the spectral approach and it is equal to the number of collocation points N coll . The expectation is that a (n) n → 1 and a (n) m → 0 if m ̸ = n. As we can see from Figs. 5 and 6, as the ratio b/a goes to zero, the convergence a (n) m → δ mn is rather good, and qualitatively improves as we refine the set of basis functions in the pseudospectral method.
In a nutshell, we see that everything in the numerical Euclidean approach is consistent with our previous real time analysis in Section 4.3.1, meaning that the quadratic kernel ∆ c diverges as L −3 when we take L → 0. As such, we conclude that the limit L → 0 of this SYM Wilson loop does not describe the physics we wish to capture, because it is dominated by UV contributions that are not in the domain of any low-energy effective description. This is also consistent with our previous discussion that we in fact expect ⟨T W [C]⟩ T = 1 for the SU(3) Wilson line configuration that is relevant to quarkonium. As such, we conclude that we must seek other configurations to describe the Wilson loop that is relevant for quarkonium dynamics in a thermal medium.

The Wilson loop with antipodal S coordinates
Having studied the supersymmetric Wilson loop with constantn, we now proceed to investigate the other candidate configuration, where we haven =n 0 on one of the timelike segments comprising the loop, andn = −n 0 on the other antiparallel segment of the loop.
This configuration has received less attention in the AdS/CFT studies of heavy quarks, mainly because it does not generate a Coulomb interaction potential between a heavy quark pair. However, while it has been usually less emphasized, it has been discussed in many AdS/CFT studies, starting from the same works that discussed the heavy quark interaction potential [86,87]. As we will also verify momentarily, it has the crucial property that ⟨W BPS [C]⟩ = 1 for a contour going from one point to another and coming back to the starting point along the same path. Apart from the fact that one can verify this identity by hand in the CFT, this relation is protected by supersymmetry [95]. Moreover, the fact that this configuration might be relevant for the dynamics of a quark-antiquark pair was hinted in Ref. [87], where the first appearance of quark pairs and heavy quark pairs featured antipodal positions on the S 5 .
Given the relative lack of attention that this configuration has received, especially for phenomenological applications, we will try to make our discussions as detailed as possible.

Background
As before, we consider the T → ∞ limit of the contour C that defines the Wilson loop from which we can extract the correlator relevant for quarkonium transport. In N = 4 SYM, we also have to specify the position on the S 5 that the Wilson loop goes over for it to have a dual gravitational description in terms of an extremal surface in AdS 5 × S 5 . As discussed in Section 3.2.1, a natural choice is to have two long timelike Wilson lines with antipodal positions on the S 5 . Without loss of generality, we can describe the distance between the S 5 coordinates on the two boundary segments by a large circle angular coordinate ϕ ∈ [0, 2π), and then we may substitute dΩ 2 5 into the metric shown in Eq. (3.16) by dϕ 2 . First we study whether there is an extremal surface connecting the two boundary segments that can be described by X µ (τ, ϕ) = (τ, 0, 0, 0, z(ϕ),n(ϕ)) , (5.1) with which the Nambu-Goto action reads: where now we define z ′ = ∂z/∂ϕ. The action of the resulting extremal surface should be compared with the action of two disconnected worldsheets falling into the black hole. If we find a positive 18 regularized action by extremizing the action (5.2), then it will be the preferred configuration, as it will be the one of the lowest energy. On the other hand, if the energy of the configuration we find by extremizing the action (5.2) is higher than  Figure 7. Left: the relation ∆ϕ(πT z m ) that determines the angular distance spanned on the S 5 by the connected configuration that reaches a maximal AdS radial coordinate z = z m . Right: dimensionless configuration energy E for the extremal worldsheet that is described by a connected configuration z = z(ϕ). ∆ϕ = ∆ϕ(πT z m ) is one to one, and that at any ∆ϕ > 0, their energy is strictly greater than that of two disconnected, radially infalling worldsheets each hanging from their respective boundary toward the bulk of AdS 5 . Crucially, this means that for our present purpose, which is to find the minimal energy configuration for ∆ϕ = π/2, the relevant extremal surface is that of two disconnected (at least at times |t| ≪ T /2) worldsheets hanging radially into the bulk of AdS 5 . The alternative solution, which is energetically disfavored, is a surface that lies at z = 0 and can be parametrized by time t ∈ (−T /2, T /2) and the angle ϕ ∈ (0, π). It is interesting to note that in the strict limit T = 0, all connected configurations have the same energy as the radially infalling solution (i.e., zero), for any value of z m . However, since we are interested in the physics in the presence of a thermal plasma, there is no ambiguity in terms of which solution to choose.
Therefore, the solution presently relevant to our purpose is parametrized by two disjoint surfaces X µ L (τ, z) = (τ, 0, 0, 0, z, −n 0 ) , (5.7) and X µ R (τ, z) = (τ, 0, 0, 0, z,n 0 ) , where z ∈ (0, (πT ) −1 ) is an independent coordinate in this description, and plays the role of one of the worldsheet coordinates. We note that the solution we found above also applies if the two timelike sides of the contour C are at nonzero spatial separation L > 0 (provided they remain at antipodal positions ∆ϕ = π/2 on the S 5 ), because allowing for a ϕ-dependent y 1 coordinate in the connected background solution can only increase the configuration energy. In summary, this configuration features two worldsheets that fall into the black hole, which intuitively represents the propagation of two unbound heavy quarks, with their interactions being screened by the thermal medium [56,62,113]. By construction, this configuration has E d = 0 (after subtracting the energy associated to the heavy quark masses), as required to satisfy ⟨W BPS [C]⟩ = 1.

Fluctuations
To calculate the response kernel of interest for this configuration ∆ d ij , we study the dynamics of fluctuations on top of the background worldsheet we just found. In consistency with the preceding discussion, we work in the limit T → ∞ (concretely, |t 1 |, |t 2 | ≪ T ), which also simplifies the calculations because of the time translational invariance. In this setup, whether T is finite or infinite does not affect the final result, as long as the time domain of interest is covered by the timelike Wilson lines, since the parts of the timelike Wilson lines that are out of the time domain of interest cancel due to U U −1 = 1. 19 As such, taking T → ∞ is not an approximation, but rather, it is manifestly equal to the starting point.
As discussed in Section 3.1, the appropriate boundary deformations of the contour C to evaluate the chromoelectric field correlator in the pure SU(N c ) gauge theory are the antisymmetric ones, as shown in Eq. (3.8). By the same argument as in that Section, these are also the only nontrivial deformations from which we can extract the desired correlation function in our present setup. The reason is that symmetric boundary deformations of the contour C preserve the value of the (supersymmetric) Wilson loop W BPS [C] = 1, which is also a consequence of having each antiparallel timelike Wilson line with antipodal positions on the S 5 .
Having made these remarks, we now proceed to introduce perturbations on top of the background worldsheet to enable us to evaluate the path functional derivatives on the boundary. Compared to the setup in the previous section, the parametrization of the perturbed worldsheet here is remarkably simpler: X µ R (t, z) = (t, +y 1 (t, z), +y 2 (t, z), +y 3 (t, z), z, +n 0 ) , (5.10) where we have already set to zero the fluctuations corresponding to reparametrization invariance, y 0 and y 4 . We have also already incorporated the fact that, because of the antisymmetry of the boundary conditions we will use, the solutions for the fluctuations on either surface will be equal but with opposite signs. The action, up to quadratic order for the fluctuations, takes the form where we have subtracted the action corresponding to the energy of two heavy quarks at rest, which is incidentally equal to the background action in this case. The form of the action for the fluctuations is the same for all components and is given by where a prime denotes the derivative with respect to z and a dot represents the derivative with respect to t.
Integrating by parts, and using the equations of motion, we can evaluate the on-shell action as which involves two explicitly nontrivial limits. Let us first focus on the first limit, where z → (πT ) −1 . While it is tempting to conclude that it is zero because f (z = (πT ) −1 ) = 0, we actually have to solve for the mode functions first and verify that its product with y ′ y indeed goes to zero. In the following Section 5.3, after selecting appropriate boundary conditions at the black hole horizon, we will confirm that this is the case. Therefore, we shall drop this term in the remainder of this section. The second limit, where z → 0, is even more subtle, because it can be manifestly divergent if y ′ does not go to zero fast enough. Nonetheless, after solving for the mode functions and investigating them at small z, we shall see that it contains a 1/z divergent term of the form discussed around Eq. (3.3) and in the footnote. 2 That is to say, it is generated by the contact term that comes from evaluating the variational derivatives with respect to the path deformations f µ at the same point. As we will see later, this contribution has exactly the same form and value as that calculated in the previous Section 4. 3. As such, we justify it to simply subtract the second limiting term from our final result. Furthermore, using that we will find the relevant mode functions satisfy y ′ (t, z = 0) = 0, we can conclude by repeated use of the L'Hopital's rule that As such, the response function we will need to calculate has a third derivative. The other way of distributing the three derivatives gives terms of the form ∂ 2 y ∂z 2 ∂y ∂z , which vanish because the mode functions satisfy y ′ (t, z = 0) = 0.
Let us now define the response function we will calculate. As in the case of the calculation of Section 4.3, we identify the value of the fluctuation on the boundary with that of the contour deformation y(t, z = 0) = h ⊥ (t). As such, we introduce the response kernel with which we find With this, the correlation function we seek is determined by As such, all that remains to be done is to evaluate the response function K d ⊥ .

Calculation of the time-ordered non-Abelian electric field correlator
To calculate the response function K d ⊥ , we proceed by varying the action S (2),⊥ NG,d with respect to y to obtain its equations of motion, and then transform to the frequency domain. Then, introducing ξ = πT z, the equation we want to solve is which is actually equivalent to the one found by Ref. [65] to calculate the heavy quark diffusion coefficient in strongly coupled N = 4 SYM theory. For the benefit of the reader, we note that in their notation, the independent variable that parametrizes the worldsheet is u = ξ 2 .
To find the solutions to Eq. (5.18), we proceed as in Ref. [65] to factor out the highly oscillatory piece that is generated close to the black hole event horizon. This can be done by the same WKB-type analysis we carried out in the calculation of the response functions in Section 4.3. Using the same notation, we introduce where the prefactor (1 − ξ 4 ) ±iΩ/4 is obtained by direct integration of the WKB phase in Eqs. (4.38) and (4.39). To facilitate comparison, we have written y ± ω in the way of Eq. (5.19) such that the resulting mode functions are the same as in Ref. [65]. With this definition, The purpose of extracting the highly oscillatory phase from y ± ω (ξ) is to get an equation for F ± ω (ξ) such that a regular solution can be found at ξ = 1. This condition must be imposed by hand, because Eq. (5.20) has two independent solutions: one regular at the horizon, and the other oscillating twice as fast as the solutions for y ± ω (ξ). Examining the differential equation for F ± ω and demanding regularity at the horizon, we find this implies where the last condition fully determines the mode solution, up to an overall normalization. This condition allows one to find numerical solutions to Eq. (5.20) ensuring regularity at the horizon. The other input required to determine the correlation function is the boundary conditions, i.e., the prescription to select the appropriate linear combination of the mode functions that determines the response kernel. Because we have extended our contour to infinity by taking the limit T → ∞, the boundary conditions are determined by the timeordering prescription ω → ω(1 + iϵ), which is a consequence of the aspects discussed in Section 3.3.2. For concreteness, let us focus on the case ω > 0. This is without loss of generality, because we are calculating a time-ordered correlation, and so, the full result will be immediately obtained by taking ω → |ω|.
With these preliminaries, we can now analyze the mode functions (5.19) under an infinitesimal complex rotation ω → ω(1 + iϵ). To select or discard a solution, we need to know whether one of the modes generates a divergent limit in the action (5.13). In particular, whether the limit exists for the mode solutions y ± ω when the iϵ prescription is taken into account. The reason why this particular limit is relevant is because of its explicit appearance in the expression for the on-shell action (5.13), which must be finite for the action to be at a well-defined extremum.
As discussed before, F ± ω is regular and finite at the horizon ξ = 1, and by inspecting the differential equation (5.20) that defines it, it is also analytic in ω. As such, no singularity will appear in F ± ω (ξ = 1) by rotating the frequency ω by a small amount from the real axis into the complex plane. Therefore, the deformation by iϵ affects the result predominantly through the WKB factor exp(± iΩ−Ωϵ 4 ln(1 − ξ 4 )). However, this means that y + ω will grow as e ϵΩ| ln(1−ξ 4 )|/4 close to the horizon for ω > 0, and therefore, substituting the mode function y + ω into Eq. (5.22) leads to a divergent limit. We then conclude that we must keep only y − ω as the allowed mode solution for ω > 0. By extension, the mode solution to be kept at arbitrary ω is y − |ω| . Now we can evaluate the response function K d ⊥ (ω) by means of substituting y − ω into Eq. (5.15). Given what we just showed, the WKB factor of y − ω goes like exp(− iΩ−Ωϵ 4 ln(1 − ξ 4 )), which goes to zero for ω > 0 as ξ → 1. Consequently, the first term of the on-shell action (5.13) vanishes, and we are left with the second term only. By direct inspection of the mode equation (5.20), we see that regularity of F ± ω requires which verifies our earlier claim that y ′ (t, z = 0) = 0. The other claim of the previous section that we have yet to verify is that the (divergent) contact terms are the same as in the heavy quark interaction potential case. To show this, we may write the unregularized response kernel K d 0⊥ from the second term in Eq. (5.13), and find where we have used L'Hopital's rule to obtain the second equality. By virtue of the second condition in Eq. (5.23), ∂ 2 y − ω ∂z 2 (z = 0) = ω 2 y − ω (z = 0), we can add and subtract the divergent part to get the unregularized response kernel as By comparison with Eq. (4.58), and keeping in mind that the z 2 m / √ f m in that expression is cancelled by the relative prefactor in the definition of ∆ c ij , it is clear that the nature of the last term in Eq. (5.25) is that of a contact term.
Therefore, switching back to ξ = πT z, we may write the regularized response kernel (which is the one that enters the chromoelectric field correlator) as The last equality follows from direct inspection of the mode functions shown in Eq. (5.19), and the fact that the prefactor (1 − ξ 4 ) −iΩ/4 has no effect in the result because all of its first three derivatives with respect to ξ vanish. In terms of ∆ d ij , the result is This result may be evaluated numerically after plugging the boundary condition shown in Eq. (5.21) to the differential equation (5.20), which defines F − ω . 20 We plot the result for a general temperature in the next section, where we give the final result of our calculation.
Before proceeding further, it is also instructive to evaluate the zero temperature limit of our expression. To do so, it is most convenient to go back to the original AdS coordinate z in the mode equation for y ω . Then, setting T = 0, Eq. (5.18) becomes ∂ 2 y ω ∂z 2 − 2 z ∂y ω ∂z + ω 2 y ω = 0 (5.28) and the solutions are relatively simpler: We choose the sign labelling so that the solutions match their finite-temperature counterparts in the limit T → 0. For visual clarity, we show the mode solutions as a function of z for a range of frequencies rescaled by the temperature in Fig. 8. The fact that we have an explicit expression then allows us to evaluate Eq. (5.26) explicitly, obtaining  We note that the cubic power in the frequency is exactly what one expects from dimensional analysis of the correlation function we are interested in from the field theory perspective. In terms of ∆ d ij , we have In summary, we have obtained the response kernel K d ⊥ that determines the on-shell Nambu-Goto action up to second order in the contour deformations of the Wilson loop expectation value that is dual to it. Contrary to the results we found in Section 4, these are well-behaved and so provide a quantitative description of the dynamics of in-medium quarkonium. From our discussion in Section 3.2.1, the background configuration is also well-founded. Therefore, we conclude that this is the N = 4 observable that most closely resembles the analogous QCD correlation function, and will use it as the N = 4 result for the quarkonium transport coefficients. 21 We give the expressions and plots for the chromoelectric field correlator that we extract from this holographic calculation in the next section, where we also discuss its implications as a baseline for phenomenological applications.  Figure 9. Real part (left) and imaginary part (right) of the non-Abelian electric field correlator of interest. The finite temperature result is shown in solid lines, and the zero temperature limit is shown in black dashed lines. The leading low-frequency limit is shown in red dotted lines. As before, the arguments of the functions at zero temperature have been rescaled by πT to have a clean visual comparison.
6 Results and applications

Main result
Since we have calculated ∆ c,d ij (ω, L) for each configuration, we are ready to conclude and give the expression for the analog object in N = 4 SYM to the time-ordered correlator of chromoelectric fields dressed by Wilson lines in QCD.
Inspecting the behavior of ∆ c ij (ω, L) as L → 0, we conclude we must discard the configuration with constantn on the two timelike segments of the Wilson line on the AdS boundary, as it does not have a sensible limit and does not satisfy ⟨W BPS [C]⟩ = 1, which is required to give the Wilson line configuration between the two non-Abelian electric fields the interpretation of an adjoint Wilson line. On the other hand, the configuration where the two timelike Wilson lines have antipodal positions on the S 5 does provide a sensible answer with perturbations that modify the energy of the worldsheet by a finite amount, and moreover, it fulfils all of the expectations that we require for the standard QCD Wilson loop. This is further substantiated by our discussion in Subsection 3.2.1, where we advanced that the appropriate description for quarkonium transport in medium is given by the Wilson loop wheren takes antipodal positions on the S 5 for the two timelike segments. As such, we use the latter one and obtain: This is the main result of this paper. We plot this in Fig. 9. Furthermore, its zerotemperature limit is given by The other limit of interest is the low-frequency limit, which can be extracted analytically by solving the mode equation (5.20) up to linear order in Ω. The algebraic steps necessary to do this are the same as those in the heavy quark diffusion coefficient calculation [65]. The result is where we have kept one higher order in ω/T than that explicitly shown in Ref. [65]. The details of how this expansion was carried out can be found in Appendix D.1. With the expression for the correlation function in hand, we can now use it to describe how a heavy quark-antiquark pair will propagate through the thermal N = 4 SYM plasma. Specifically, we can calculate the spectral function that determines the transition rates of in-medium quarkonium within a potential non-relativistic EFT description as shown in Section 2, and draw the phenomenological implications thereof for strongly coupled plasmas. This is what we will first show in the next section. As an additional application of our result, we will compare our result to its weak-coupling limit in N = 4 SYM, and lay out the case for computing subleading corrections, both at weak and strong coupling.

Applications
Having calculated the non-Abelian electric correlator of interest, we can now discuss the physical consequences and applications of our result. We shall focus on two pieces of theoretical understanding we can gain from here. First, we will discuss the immediate consequences of our result for in-medium quarkonium dynamics in a strongly coupled N = 4 supersymmetric Yang-Mills plasma. After that, we will discuss how a comparison with the weak coupling expansion of the same correlator may be used together with our strong coupling result to advance our understanding of the always elusive non-perturbative dynamics of gauge theories at intermediate coupling.

Evaluation of chromoelectric field spectral function
We are now ready to evaluate the spectral function that encodes important information of the plasma relevant for quarkonium transport, by using the relations introduced in Section 2, as dictated by our N = 4 SYM result. Specifically, we use In general, the above equation would be an integral expression that can only be evaluated numerically. However, in N = 4 Yang-Mills theory, the fact that one obtains [g T E ] N =4 (ω) by selecting modes using the time-ordering prescription ω → ω(1 + iϵ) gives us more analytic control. As we show in Appendix D.2, this property allows us to prove that and consequently, the correlation function that enters the quantum and classical quarkonium time evolution equations is given by  Figure 10. Spectral function for quarkonium transport. Only the positive frequency domain is shown, as ρ ++ E vanishes for ω < 0.
One more step gives us an explicit expression for the spectral function which we plot in Fig. 10. This function is manifestly neither even nor odd, as expected from the evidence coming from the perturbative calculations [45,47]. One immediate implication of our results, which may already be seen from Fig. 9 is that the transport coefficients introduced in the Quantum Brownian motion limit [23,27,34,35,46] of the open quantum system approach to in-medium quarkonium, namely, the analogs to vanish in strongly coupled N = 4 supersymmetric Yang-Mills theory. Explicitly, The quantity γ adj represents the mass shift of quarkonium states inside a plasma. The result γ N =4 adj = 0 is consistent with a recent lattice QCD study [76]. In the quantum optical limit where the quarkonium time evolution can be effectively described by a Boltzmann equation as shown in Eq. (2.10), it is the finite frequency part of the chromoelectric field correlator that enters the quarkonium dissociation and recombination rates. Because the argument of [g ++ E ] > is negative in Eq. (2.11), our result in Eq. (6.6) indicates that the dissociation rate of a small-size quarkonium state in a strongly coupled QGP vanishes. Using Eq. (2.21), we also see that the recombination rate in Eq. (2.12) also vanishes in the same limit.

Towards intermediate couplings in N = 4 SYM
A natural question we can ask is how the calculation result in the strong coupling limit compares with that in the weak coupling limit. We also want to understand if they allow for an interpolation at intermediate couplings.
At weak coupling, the non-Abelian electric field correlation function we have set out to calculate can be evaluated directly using the standard real-time perturbation theory in the Schwinger-Keldysh formalism. In the large N c limit of N = 4 SYM, it reads It is apparent that the small ω/T limit of the strong coupling and weak coupling results is different: in the weakly coupled case it goes as |ω| 2 , while for the strongly coupled limit it is linear in |ω|. That is to say, for the range of frequencies that is sensitive to thermal effects, the physics at weak and strong coupling is different. Note that the |ω| 2 behavior displayed above implies that the transport coefficients κ N =4 adj and γ N =4 adj vanish at leading order in perturbation theory. At NLO, however, κ N =4 adj is known to be nonzero [116], and is actually equal to κ N =4 fund up to this order in perturbation theory. Therefore, it seems that the most interesting thermal physics lies in the intermediate coupling regime. However, a first approximation to this regime via interpolation between weakly and strongly coupled results would require to calculate the first nonvanishing contributions from both sides, which is a challenging computation we do not undertake in this work.
On the other hand, in the T = 0 limit, the frequency dependence of both results agrees: both are proportional to |ω| 3 . This is unsurprising given that N = 4 supersymmetric Yang-Mills theory is a conformal field theory, and there is thus no other scale available to give rise to a different behavior. As such, in vacuum we have We plot both limits in Fig. 11, together with the Padé approximant of order [2/1] in √ λ that interpolates between the two limits. Such an interpolation constitutes, at most, an educated guess of the result for the chromoelectric correlation function in the intermediate coupling regime around λ ∼ O (10). As it should be clear from the comparison, we only expect the result to be valid asymptotically.
Nonetheless, such a comparison may provide valuable insight into what the behavior of the correlator is at intermediate couplings. In fact, Fig. 11   is a next-leading order calculation analogous to what has already been done for QCD in Ref. [45], but this time for N = 4 SYM. At strong coupling, one would have to evaluate the quantum corrections to the string worldsheet action in the so-called semiclassical expansion, which is tantamount to a 1-loop calculation of the fluctuation fields on the worldsheet [100]. Both are necessary steps towards a more complete understanding of the correlator, which are along the path that we want to follow in the future (in the hope that the convergence of the series is comparable to that of other observables in N = 4 SYM, e.g., the thermodynamic pressure [117]).

Conclusions
In this paper we studied the gauge-invariant correlation function of two chromoelectric fields connected via an adjoint Wilson line at finite temperature. In QCD, this correlation function determines the in-medium dynamics of small-size quarkonium states. However, in light of the tools provided by the gauge/gravity duality at strong coupling, and of the impressive success that holography has had in describing the properties of the quark-gluon plasma, we studied this correlation function in N = 4 supersymmetric Yang-Mills theory. By expressing the correlation of two non-Abelian electric fields at different times connected via Wilson lines as a variation of a fundamental Wilson loop, we were able to calculate it by evaluating the variations of the holographic dual description of the Wilson loop, which is given by an extremal surface in AdS 5 × S 5 .
We examined two candidates for the N = 4 Yang-Mills description of the relevant Wilson loop using the AdS/CFT correspondence. We considered two worldsheet configu-rations: one with a constant S 5 coordinate on the two timelike segments of the Wilson loop on the AdS boundary and the other with antipodal S 5 coordinates. In our calculations, we first obtained the background worldsheet for each configuration, which is an extremal surface that hangs from the Wilson loop contour on the AdS boundary. Then we introduced deformations of the Wilson loop contour and studied how the perturbations propagate on the worldsheet. Finally we obtained the correlator of our interest by calculating the derivative response function on one side of the timelike segments of the Wilson loop caused by the deformation on the other.
We showed that the non-Abelian electric field correlator obtained from the configuration with a constant S 5 coordinate diverges and does not give a sensible result for the correlator of our interest. This is related to the fact that this configuration does not satisfy our expectation about the pure gauge Wilson loop ⟨W [C]⟩ = 1, where C denotes a Wilson loop consisting of two antiparallel timelike Wilson lines on top of each other. On the other hand, the configuration with antipodal S 5 coordinates is well-motivated from a dynamical perspective, as it describes the propagation of two heavy probes with conjugate (opposite) charges. Furthermore, it is consistent with the properties required for a dual description of the pure gauge Wilson loop according to the prescriptions that have been proposed for this type of loop, and leads to a finite and sensible result for the correlator. Following these observations, we conjectured that the pure gauge adjoint Wilson line in the strong coupling limit of N = 4 SYM may be described by two antiparallel 1/2 BPS Wilson lines with antipodal values of the S 5 coordinate. Further verification of this proposal, such as evaluating this configuration for a finite temporal extent (as opposed to the limit we considered in this paper, where we let T → ∞), is an interesting endeavor that we hope to undertake in the near future.
Our results indicate that the two transport coefficients in the Lindblad equation for quarkonium in-medium dynamics in the quantum Brownian motion limit, as defined through the chromoelectric field correlator we studied, vanish in a strongly coupled N = 4 SYM plasma. Furthermore, the correlation functions [g ±± E ] > that determine quarkonium dissociation and recombination in the Boltzmann equation that is valid in the quantum optical limit also vanishes in a strongly coupled N = 4 SYM plasma. This means the in-medium dynamics of small-size quarkonium states is trivial at leading order in the multipole expansion, which is used to simplify the evolution equations in both the quantum Brownian motion and quantum optical limits. The reason behind this triviality can be multi-fold. Firstly, in the strongly coupled, N c → ∞ limit of N = 4 SYM theory, the energy gap between a heavy quark-antiquark pair in the color singlet and an octet pair, as extracted from the expectation value of a rectangular Wilson loop, is infinitely large. This means that it would take an infinite amount of energy to break the bound color singlet pair, which a QGP at finite temperature does not have. The inverse process, i.e., recombination, cannot happen either, because the process would release an infinite amount of energy, which cannot be absorbed by the local degrees of freedom of a QGP at finite constant temperature. Secondly, the validity of both the quantum Brownian motion and quantum optical limits may break down when the plasma is strongly coupled, which should be further investigated in the future. Given the fact that some properties of the QGP at low temperatures can be well described by a strongly coupled N = 4 SYM plasma, going beyond these limits might prove of paramount relevance for phenomenology.
Furthermore, it is unclear to what extent the vanishing result may be a consequence of the leading description in √ λ. This motivates going beyond the leading strong coupling limit in N = 4 SYM by accounting for effects of quantum fluctuations on the worldsheet on the gravitational side of the duality. It also motivates calculating the subleading terms in a weak-coupling perturbative expansion in order to have an improved Padé approximation of the complete correlation function, both for the large frequency regime as well as for the low frequency regime. Once calculated, these results would provide insight into the intermediate coupling regime, and therefore improve our description of quarkonium dynamics by setting the coupling at phenomenologically relevant values.
On the QCD side, our results further advance the need for a non-perturbative determination of the chromoelectric correlator that determines in-medium quarkonium dynamics. Current and past estimates of the transport coefficients rely on spectral function reconstructions determined either from the correlator that describes heavy quark diffusion or from the heavy quark interaction potential dynamics. Neither provides direct access to the correlation function we have discussed in this paper. Given that the results for the heavy quark diffusion coefficient κ and the quarkonium transport coefficient κ adj differ in N = 4 SYM, the expectation is that this will also be the case in QCD. Therefore, a non-perturbative calculation of the quarkonium transport coefficients has the potential to make a significant contribution to phenomenological studies of in-medium quarkonium.
Finally, from a data-driven perspective, one should try to extract this correlation from phenomenological studies and experimental measurements by applying Bayesian analysis techniques, where one uses some ansatz for the correlation that is well-motivated from both weak-coupling and strong-coupling studies, varies the parameters in the ansatz and compares the calculation results of quarkonium nuclear modification factors and elliptic flow coefficients with experimental data. The Bayesian analysis is then applied to systematically find the best set of parameters in describing the data and estimate the parameters' uncertainties. For example, it will be particularly instructive to see how data from 200 GeV Au-Au collisions, particularly from the sPHENIX program, can be used to better constrain the finite frequency dependence of the correlator via a Bayesian analysis. We expect data coming from these collision energies will be sensitive to its finite frequency dependence, because the prevailing temperature regime in this experiment is of low temperature comparable to the energy level splittings of quarkonia states. All these studies will deepen our understanding of quarkonium in-medium dynamics and the relevant transport properties of the QGP.

A KMS relations of chromoelectric correlators
In this Appendix we will verify the KMS relation between [g ++ E ] > (t) and [g ++ E ] < (t), as introduced in the main text in Eqs. (2.14) and (2.15). For the proof of how the timereversal symmetry relates [g ++ E ] and [g −− E ], we refer the reader to Appendix A of [45]. The following proof adds an omission that two of us made in Appendix A of [45], where we did not pay enough attention to the effects of the adjoint color charge on the thermal average of the medium. We hope that the subsequent discussion makes these aspects more transparent.
We begin by noting that an adjoint Wilson line in the interaction picture can be written in terms of time-evolution operators as where H is the QGP Hamiltonian, and one may interpret Hδ ab − gA c 0 (0)[T c Adj ] ab as the total Hamiltonian when there is a point color charge in the adjoint representation at the position x = 0.
With this, the Wightman correlator [g ++ E ] > (t) can be written as The KMS conjugate of this correlator is obtained by shifting t → t − iβ. We therefore obtain We can explicitly see in this expression that the thermal ensemble is now determined by the total Hamiltonian Hδ ab − gA c 0 (0)[T c Adj ] ab instead of just H. The equivalence with Eq. (2.15) is then verified by using that

B.1 Time-ordered products of Wilson lines
In this work, we deal with the calculation of the time-ordered correlator where W ab [t,0] is an adjoint Wilson line, that is written in terms of fundamental Wilson lines as where the time-ordering symbol is necessary to preserve the explicit ordering of operators in an adjoint Wilson line. For concreteness, we write both adjoint and fundamental lines below: where the difference is the representation of the SU(N c ) generator matrices. Note that the operator ordering in the adjoint Wilson line is not the "natural" one in terms of fundamental Wilson lines, because the operator products are not ordered in the same way as the matrix products that contract color indices. In this sense, Eq. (B.2) without the time-ordering symbol is only indicative of the color product structure, but is not an explicit expression in terms of how the gauge field operators A a 0 (t) are ordered. This is specified by the symbolT , but by itself does not provide an explicit method to calculate it. One explicit way to evaluate it is given by its path integral representation, to which we will explain in a moment. Before doing that, however, it is useful to discuss how this time ordering of fundamental Wilson lines appears from the dynamics of two (coincident) point color charges, which we take to be in the fundamental and anti-fundamental representations of SU(N c ).
To accomplish this, let us collectively denote the "colors" of the QQ pair by QQ ij , with i being the index of the quark in the fundamental representation, and j the index of the anti-quark in the anti-fundamental representation. The dynamics are given by where we have used that T a Anti−Fund ij = −[T a Fund ] ji .
Formally, we can write the solution to this equation as QQ ij (t) = W ii 0 ,j 0 j (t) QQ i 0 j 0 (t = 0) , (B.6) where W ii 0 ,j 0 j (t) obeys the same equation as QQ ij with W ii 0 ,j 0 j (t = 0) = δ ii 0 δ jj 0 as the initial condition. Note that, by construction, we have A quick way to see this is to note that if A a 0 were ordinary numbers, then the time ordering would be irrelevant and the Wilson lines in the last expression would be completely decoupled from each other. However, because A a 0 (t) are in principle non-commuting operators, we have to keep track of the fact that A a 0 (t) is always inserted to the left of operators A a 0 (t ′ ) when t > t ′ . This is due to A a 0 (t) being to the left of W in its defining differential equation (B.7). Therefore, what we get out of W ii 0 ,j 0 j is a fundamental and an anti-fundamental Wilson line put together, with their operators time-ordered.
A consistency check is to verify that we can get the adjoint Wilson line from W ii 0 ,j 0 j . Indeed, we can consider and get, after contracting the color indices and using the Lie Algebra of the group, which is exactly the defining equation for an adjoint Wilson line: It is also direct to see from here that the Wilson loop we consider in the main text satisfies and contracting the indices in the previous equations, we see that and therefore, given the initial condition W ii 0 ,j 0 j (t = 0) = δ ii 0 δ jj 0 , we conclude that as claimed in the main text. It is worth noting that this is self-evident from the path integral formulation. Indeed, collectively denoting the field content of the theory by φ, one can calculate time-ordered (vacuum) correlation functions as The step to the last line is achieved because, inside the path integral, the Wilson lines are just SU(N c ) unitary matrices that are inverses of each other. The fact that we get one out of this is evidently consistent with our previous discussion.

B.2 Standard products of Wilson lines
In other contexts, it is also possible for the Wilson  which also equals unity. The reason for that here, however, is that the operators U [t,0] and U † [t,0] are inverses of each other, and should be interpreted as written, with the operator products appearing in the same way as the color products.
Interestingly, the path integral description of this object is less simple than for the time-ordered loop. The reason for this is that inserting complete bases of states along the operator products to convert the expectation value into a path integral requires following the time contour defined by the explicit operator ordering in the correlation function. In the time-ordered case, operators are, by definition, arranged further to the left at later times, and hence it is sufficient to insert complete bases of states that span the [0, t] time interval once. However, for the loop considered in this section, one has to insert complete bases of states along the interval [0, t] once in the forward direction (for U ), and once in the backward direction (for U † ).
In the context of the AdS/CFT correspondence, this type of operator ordering is realized by using the gauge/gravity duality for each segment of the path integral contour and imposing appropriate matching conditions [93,94,104]. Both the heavy quark diffusion coefficient calculation [65] and the jet quenching parameter claculation [64] implicitly have this feature. This is in contrast to the calculation presented herein, which does not require to match the background solution across manifolds that have different segments of the Schwinger-Keldysh contour as their boundaries.
C Calculation of the transverse kernel with constant S 5 coordinate In this section we evaluate the correlation kernels K AB,c ⊥ . The first step is to take the limit ξ → 0, as this is part of the definition of the correlator at any L. We use the notation F ± ω with the understanding that its dependence on ω will often be through Ω = ω/(πT ). We start from the expressions (4.57) and (4.58).
The relevant limit corresponds to taking z m → 0, which is equivalent to taking L → 0 at fixed T . Furthermore, because we have an explicit factor of z −3 m , we have to calculate the series expansion of the rest of the expression up to order z 3 m , so that in the end we get a result of the form Let us examine this term by term. We first note that we need to evaluate all terms in this expression up to cubic power in z m . In particular, we have to evaluate at least up to cubic order in z m . As with the rest of this expression, this requires to solve for F ± ω up to O(h 3 ), where h = πT z m . However, by simple inspection, one quickly realizes that the structure of the solution, up to O(h 3 ), is of the form , (C.4) which implies that the product F − ω (ξ)F + ω (ξ) is an even function of h (similarly for F + ω +F − ω ). This in turn implies that the whole object is an odd function of h. Moreover, since no terms in the power series up to O(h 3 ) involve the temperature explicitly, we have that both correlators are of the form In terms of the time coordinate, both terms only contribute to the infinitesimal neighborhood of t − t ′ = 0, and therefore we anyways expect a divergence. The leading nontrivial dependence on ω/T appears at linear order in L, which is to say, from the O(h 4 ) corrections from each term. For future reference, note that F (0) ω , and F (2) ω can be determined explicitly: where ϕ (1) ω = 2 ω (ξ) and ϕ (5) ω have a nontrivial dependence on ω. We can absorb the iϵ into the definition of ω and rotate ω → ω(1 + iϵ) at the end. It means we can expand the trigonometric functions in the definitions of the correlators and proceed without obstacle. However, because the structures are slightly different, we proceed separately for each correlator.
and instead of attempting to find a solution for arbitrary Ω, let us expand the solution in powers of Ω = ω πT . To that end, we write Higher order terms may be obtained by solving the differential equation (D.1) up to higher powers of Ω.

D.2 Consequences of the pole positions of the time-ordered correlator
In this section we shall prove that [g T E ] satisfies Eq. (6.5). Up to overall factors, and setting the normalization F ± ω (ξ = 0) = 1, the time-ordered correlator we obtained is given by (D.8) which we obtained by shifting ω → ω(1 + iϵ), which is essentially a prescription to avoid potential poles along the real ω axis.
In order to prove that [g T E ] satisfies Eq. (6.5), we will take a seemingly disconnected starting point, which nonetheless will allow us to prove our claim. To begin, consider the integral I(ω) = ∞ 0 dp 0 2ωG(p 0 ) p 2 0 − ω 2 + iϵ . (D.9) Specifically, we will prove that Im {I(ω)} = 0. Note that this integral only involves ω > 0. Then, because we constructed G(ω) by shifting potential poles on the positive real axis towards the lower half of the complex plane, there is no obstruction to Wick-rotate the integral onto the positive imaginary axis. This is possible because F − ω itself is an analytic function, provided we handle the potential UV divergences properly. We will deal with the potential large ω divergences in the next subsection.
After doing the Wick rotation, we get We have also verified this relation numerically for the vacuum-subtracted correlation functions (i.e., for ∆G(ω) = G(ω) − G(ω) T =0 ), where all integrals are convergent. For the vacuum part, where the integrals are UV-divergent because of the ω 3 power-law behavior, we present a proof in the next subsection.

D.3 UV divergent pieces in the Wick rotation
To be sure that we have the correct expression for all contributions in the above, we may work out the contributions proportional to ω 3 in the correlator independently. Because the integrals are divergent, and we do not have a natural regulator that respects all of the AdS symmetries, we will use Lorentz covariance of the boundary theory and the fact that N = 4 SYM is a conformal field theory (CFT). Once firmly on the side of the boundary theory, we may use all of the standard dimensional regularization machinery to calculate the integrals.
At T = 0 we have restored Lorentz covariance of the boundary theory, and therefore we can obtain the same correlation function but with the Wilson lines at an angle with the time axis by applying boosts. This, plus the fact that the theory is a CFT means that G(ω) ∝ |ω| 3 may be derived by integrating a momentum-space two-point function of a massless particle: (D.14) Now, we simply verify the right hand side of Eq. (D.13) mode by mode, i.e., sgn(ω)Re iω 2 ω 2 − k 2 + iϵ = ∞ −∞ dp 0 πp 0 Im i(ω + p 0 ) 2 (ω + p 0 ) 2 − k 2 + iϵ .
(D. 15) This identity is indeed satisfied, because we may write the numerator of the integrand on the right hand side as (ω + p 0 ) 2 = (ω + p 0 ) 2 − k 2 + k 2 : the first two terms then cancel the denominator, leaving their contribution as the Cauchy principal value integral of 1/p 0 , which vanishes. The last term gives a contribution that can be cast in the form of a Dirac delta function by means of ∞ −∞ dx P 1 (x − 1)(x 2 − a 2 ) = π 2 2 δ(|a| − 1) , (D. 16) and the left hand side may be immediately seen to be proportional to a Dirac delta function δ(ω 2 − k 2 ). Verifying that the coefficients match is straightforward. Therefore, the zero-temperature piece of [g T E ] satisfies Eq. (6.5), as does the thermal contribution. Thus, we have verified the claim presented in the main text.