---
abstract: 'Exact sum rules for the longitudinal and transverse part of the vector channel spectral functions at nonzero momentum are derived in the first part of the paper. The sum rules are formulated for the finite temperature spectral functions, from which the vacuum component has been subtracted, and represent a generalization of previous work in which sum rules were derived only for the zero-momentum limit. In the second part of the paper, we demonstrate how the sum rules can be used as constraints in spectral fits to lattice data at various temperatures, with the latest dynamical lattice QCD data at zero momentum.'
author:
- Philipp Gubler
- Daisuke Satow
title: Finite temperature sum rules in the vector channel at finite momentum 
---

Introduction {#sec:intro}
============

Hadrons, whose dynamics is governed by quantum chromodynamics (QCD), are deconfined at high temperature ($T$), where quarks and gluons are expected to be the fundamental degrees of freedom. Such matter is called quark-gluon plasma, and is experimentally investigated by heavy ion collision experiments. In analyzing the experimental data, electromagnetic (EM) probes such as dilepton spectra are particularly useful [@Adare:2014fwh] because once generated in the medium, they are expected to reach the detector without further QCD interaction with other particles. The electric conductivity is also an important quantity since it may increase the lifetime of magnetic fields generated in the early stage of the heavy ion collision [@McLerran:2013hla; @Tuchin:2013apa; @Tuchin:2015oka; @Hattori:2016emy]. Furthermore, how the spectrum of the vector meson is modified at finite $T$ has been discussed for a long time from the point of view of the chiral symmetry restoration [@Bernard:1988db]. The spectral function of the EM current at finite $T$, which is the focus of this paper, contains information on all the above three quantities.

There are many approaches for evaluating the EM spectral function at finite $T$, such as perturbative QCD [@Baier:1988xv], holographic QCD [@CaronHuot:2006te], model calculations [@Gale:2014dfa], sum rules [@Gubler:2015yna; @Huang:1994fs; @Kapusta:1993hq; @Zschocke:2002mn; @Kim:2017nyg], low-energy effective theory [@Chanfray; @Klingl], and lattice QCD [@Ding:2010ga; @Meyer:2011gj; @Brandt:2015aqk; @Amato:2013naa; @Aarts:2007wj]. Nevertheless, none of these approaches are perfect. Especially in the lattice QCD approach, which allows fully nonperturbative first principle QCD calculations, the spectral function cannot be analyzed directly since it is a quantity defined in real, not imaginary time. Therefore, one needs to make some assumption about the functional form of the spectral function to analyze it, or otherwise has to rely on some method to analytically continue imaginary time data to real time, such as the maximum entropy method (MEM) [@Jarrell:1996rrw; @Asakawa:2000tr; @Gubler:2010cf], the Backus-Gilbert method [@Backus:1968; @Backus:1970; @Brandt:2015sxa] or the Schlessinger point method [@Schlessinger:1967; @Tripolt:2016cya] (see Ref.[@Tripolt:2017] for a comparison of these three methods).

In our previous work of Ref.[@Gubler:2016hnf], we derived three sum rules at zero momentum, which constrain the spectral function, and used them to improve the ansatz employed in previous lattice QCD analysis. One aim of this paper is to derive similar sum rules for the small but finite spatial momentum case. At finite momentum, the EM spectral function no longer has a single independent component, but two, corresponding to the transverse and the longitudinal channels. In the longitudinal channel, a novel and robust structure, a sharp peak corresponding to the diffusion mode of the EM charge, appears.

As shown in Ref.[@Gubler:2016hnf], the sum rules can be used to improve the analysis of lattice QCD data by constraining the shape of the spectral functions and derive transport coefficients that do not appear directly in the spectral functions. In Ref.[@Gubler:2016hnf], this was demonstrated by using two sum rules (1 and 3 in this and the previous work), but the other sum rule (2) was not used. The other aim of this paper is to update the analysis such that all three sum rules can be used, and to employ the latest lattice QCD data including dynamical quarks as input for the spectral function fit.

The paper is organized as follows: In the next section, we introduce the quantities in quantum field theory that are necessary in our analysis, and explain how to derive the sum rules from the operator product expansion (OPE) and hydrodynamics [@Romatschke:2009ng]. Section \[sec:sumrule\] is devoted to the derivation of the sum rules in transverse and longitudinal channels, at small but finite spatial momentum. We also confirm that the spectral function evaluated at weak coupling and in the chiral limit satisfies these sum rules, and check to which energy region the sum rules are sensitive. We demonstrate that the sum rules can be used to improve the lattice QCD analysis for the zero-momentum case in Sec. \[sec:lattice\]. We summarize the paper and give concluding remarks in Sec. \[sec:summary\]. In the three appendixes, we evaluate the contributions to the spectral functions from the transport peak, the continuum, and the UV tail, at weak coupling and in the chiral limit.

In this paper, we recapitulate known results such as the sum rules in the transverse channel at zero momentum, Eqs.(\[eq:sumrule-1-p0\]) and (\[eq:sumrule-3-T\]). Their derivation can be found in our previous work [@Gubler:2016hnf], but we rederive them to make our paper self-contained. The recapitulation of the evaluation of the transport peak, the continuum, and the UV tails in the three appendixes is provided for the same reason.

Preliminaries
=============

In this section, we explain the method for deriving the sum rules developed in Ref. [@Romatschke:2009ng]. Only the asymptotic behaviors of the EM current correlator in the UV and IR energy regions are necessary for this purpose. We also discuss these behaviors and give their analytic form obtained from OPE and hydrodynamics in this section.

Formalism
---------

We begin by introducing quantities that will be used in the derivation of the sum rules. The retarded Green function of the EM current ($j^\mu\equiv e\sum_f q_f \overline{\psi}_f \gamma^\mu \psi_f $) is defined as $G^R_{\mu\nu}(\omega,{{\mathbf{p}}})\equiv i \int dt \int d^3{{\mathbf{x}}}e^{i\omega t-i{{\mathbf{p}}}\cdot {{\mathbf{x}}}}
\theta(t) \langle [j_\mu(t,{{\mathbf{x}}}),j_\nu(0,{{\mathbf{0}}})] \rangle$, where the average is taken over the thermal ensemble, $e$ is the coupling constant of quantum electrodynamics, $q_f$ is the charge of each quark flavor (in units of $e$), and $\psi_f$ the quark field with flavor $f$, respectively.

At finite temperature, the medium effect breaks Lorentz symmetry so that the tensor structure of the Green function has two independent components, $$\begin{aligned}
\label{eq:tensor-decomposition}
G^R_{\mu\nu}(p)&=G^T(p)P^T_{\mu\nu}(p)
+G^L(p)P^L_{\mu\nu}(p),\end{aligned}$$ where $p^\mu\equiv(\omega,{{\mathbf{p}}})$ is a shorthand notation for the energy and the spatial momentum, and $P^T_{\mu\nu}(p)\equiv g^{i}_\mu g^{j}_\nu (\delta^{ij}-{p^i p^j}/{{{\mathbf{p}}}^2})$ and $P^L_{\mu\nu}(p)\equiv P^0_{\mu\nu}(p) -P^T_{\mu\nu}(p)$ are the projection operators to the transverse and longitudinal parts with $P^0_{\mu\nu}(p)\equiv -(g_{\mu\nu}-{p_\mu p_\nu}/{p^2})$. The first (second) term in Eq. (\[eq:tensor-decomposition\]) corresponds to the transverse (longitudinal) component in three dimensions. When ${{\mathbf{p}}}$ is along the $z$-direction, they are related to the components of the Green function as $G_{T}=G^R_{11}=G^R_{22}$ and $G_L=p^2G^R_{00}/{{\mathbf{p}}}^2$.

Here we recapitulate the method of deriving sum rules at finite temperature that was developed in Ref.[@Romatschke:2009ng]. Using the residual theorem for the contour[^1] drawn in Fig. \[fig:contour\], we get $$\begin{aligned}
\label{eq:residue-theorem}
\begin{split}
&\delta G^{R}_{\mu\nu}(i\omega,{{\mathbf{p}}})-\delta G^{R}_{\mu\nu}(\infty,{{\mathbf{p}}}) \\
&~~~= \frac{1}{2\pi i}\oint_C d\omega' 
\frac{\delta G^{R}_{\mu\nu}(\omega',{{\mathbf{p}}})-\delta G^{R}_{\mu\nu}(\infty,{{\mathbf{p}}})}{\omega'-i\omega},
\end{split}\end{aligned}$$ where $\delta$ stands for the subtraction of the $T=0$ part: $\delta G^R_{\mu\nu}\equiv G^R_{\mu\nu}-G^R_{\mu\nu}|_{T=0}$. Because of this subtraction of the zero temperature part and another subtraction of $\delta G^{R}_{\mu\nu}(\infty,{{\mathbf{p}}})$ done in the expression above, all the UV divergences are regularized in all the cases we consider. Thus, the integral on the contour $C$ can be safely replaced with the integral on the real axis. Moreover, in deriving Eq. (\[eq:residue-theorem\]) we have used the property that the retarded Green function is analytic in the upper $\omega'$ plane.

Now, by taking the $\omega\rightarrow 0$ limit, Eq. (\[eq:residue-theorem\]) reduces to $$\begin{aligned}
\begin{split}
&\delta G^{R}_{\mu\nu}(0,{{\mathbf{p}}})-\delta G^{R}_{\mu\nu}(\infty,{{\mathbf{p}}}) \\
&~~~= {\text P}\frac{1}{\pi i}\int^\infty_{-\infty} d\omega' 
\frac{\delta G^{R}_{\mu\nu}(\omega',{{\mathbf{p}}})-\delta G^{R}_{\mu\nu}(\infty,{{\mathbf{p}}})}{\omega'},
\end{split}\end{aligned}$$ where we have used $1/(\omega'-i\omega)\rightarrow {\text P}(1/\omega')+i\pi\delta(\omega')$. We consider only the case of $\mu=\nu$ in this paper. Then, the real (imaginary) part of $\delta G^{R}_{\mu\nu}(\omega,{{\mathbf{p}}})$ is even (odd) in terms of $\omega$. This property enables us to simplify the equation above as $$\begin{aligned}
\label{eq:derivation}
\begin{split}
\delta G^{R}_{\mu\nu}(0,{{\mathbf{p}}})-\delta G^{R}_{\mu\nu}(\infty,{{\mathbf{p}}}) 
= \frac{2}{\pi}\int^\infty_{0} d\omega
\frac{\delta \rho_{\mu\nu}(\omega,{{\mathbf{p}}})}{\omega},
\end{split}\end{aligned}$$ where we have introduced the spectral function[^2] of the EM current, $\rho_{\mu\nu}(p)\equiv{\text{Im}}G^{R}_{\mu\nu}(p)$, and changed the label of the integration variable as $\omega'\rightarrow \omega$ for simplicity. We hence see that the asymptotic behaviors of the retarded Green function in the UV and IR regions determine the integral of the spectral function.

Let us here briefly discuss the differences between the sum rules derived here and the so-called finite energy sum rules (FESR), which are widely used in the literature (see for instance Ref.[@Ayala:2013vra]). In contrast to the procedure of this paper, one in the FESR does not subtract the zero temperature contribution as we have done in Eq.(\[eq:residue-theorem\]) and below. Instead, to avoid an ultraviolet divergence, one does not take the radius of the contour in Fig.\[fig:contour\] to infinity, but sets it to some threshold value, which however should be large enough such that the OPE is still approximately valid. Doing this, one can derive sum rules that are not exact, but practically useful, as they can constrain the spectral function below the threshold value that usually contains the most interesting physical content.

![The contour $C$, used in the integral of Eq. (\[eq:residue-theorem\]).[]{data-label="fig:contour"}](contour.eps){width="25.00000%"}

UV behavior
-----------

The behavior in the UV region can be described with the help of the OPE. At leading order in the coupling constant, the Wilson coefficients of the operators with dimension 4 read [@CaronHuot:2009ns; @Shifman:1978by] $$\begin{aligned}
\label{eq:OPE}
\begin{split}
& \delta G^{R}_{\mu\nu}(\omega,{{\mathbf{p}}}) \\
&=e^2\sum q^2_f \frac{1}{p^2} \\
&~~~\times \Bigl[\left\{2m_f \delta\left\langle\overline{\psi}_f \psi_f \right\rangle
+\frac{1}{12}\delta\left\langle \frac{\alpha_s}{\pi}G^2\right\rangle \right\}P^0_{\mu\nu}(p) \\
&~~~-2\delta\left\langle T^{\alpha\beta}_f \right\rangle A_{\mu\nu \alpha\beta}(p)
\Bigr]
+{\cal O}(\omega^{-4}), 
\end{split}\end{aligned}$$ where $A_{\mu\nu \alpha\beta}(p)\equiv g_{\mu\alpha}g_{\nu\beta}+g_{\mu\beta}g_{\nu\alpha}-2(g_{\mu\alpha}p_{\nu}p_{\beta}+
g_{\nu\alpha}p_{\mu}p_{\beta}-g_{\mu\nu}p_\alpha p_\beta)/p^2$, $G^{\mu\nu}_a\equiv \partial^\mu A^\nu_a-\partial^\nu A^\mu_a -gf_{abc} A^\mu_b A^\nu_c$ is the field strength, $G^2\equiv G^a_{\mu\nu}G^{a\mu\nu}$, $T^{\alpha\beta}_{f}\equiv i{\cal {ST}} \overline{\psi}_f\gamma^\alpha D^\beta \psi_f $ is the quark component to the traceless part of the energy-momentum tensor, $D^\mu\equiv \partial^\mu +igA^\mu_a t^a$ is the covariant derivative, $A^\mu_a$ is the gluon field, $t^a$ is the generator of the $SU({N_c})$ group in the fundamental representation, $f_{abc}$ is the structure constant of the $SU({N_c})$ group, $m_f$ is the current quark mass, $g$ is the QCD coupling constant, $\alpha_s\equiv g^2/(4\pi)$, and ${N_c}$ is the number of the colors. ${\cal {ST}}$ makes a tensor symmetric and traceless, ${\cal {ST}}O^{\alpha\beta}\equiv (O^{\alpha\beta}+O^{\beta\alpha})/2-g^{\alpha\beta}O^\mu_\mu /4$.

We decompose Eq.(\[eq:OPE\]) into transverse and longitudinal components as $$\begin{aligned}
\nonumber
\delta G_{T}(\omega,{{\mathbf{p}}}) 
&=e^2\sum q^2_f \frac{1}{p^2} \\
\nonumber
&~~~\times \Bigl[\left\{2m_f \delta\left\langle\overline{\psi}_f \psi_f \right\rangle
+\frac{1}{12}\delta\left\langle \frac{\alpha_s}{\pi}G^2\right\rangle \right\}\\
\label{eq:OPE-wo-RG-T}
&~~~+\frac{8}{3}\frac{\omega^2+{{\mathbf{p}}}^2}{p^2}\delta\left\langle T^{00}_f \right\rangle 
\Bigr]
+{\cal O}(\omega^{-4}),\\
\nonumber
\delta G^{R}_{00}(\omega,{{\mathbf{p}}}) 
&=e^2\sum q^2_f \frac{1}{p^2}\frac{{{\mathbf{p}}}^2}{p^2} \\ 
\nonumber
&~~~\times \Bigl[\left\{2m_f \delta\left\langle\overline{\psi}_f \psi_f \right\rangle
+\frac{1}{12}\delta\left\langle \frac{\alpha_s}{\pi}G^2\right\rangle \right\} \\
\label{eq:OPE-wo-RG-L}
&~~~+\frac{8}{3} \delta\left\langle T^{00}_f \right\rangle  
\Bigr]
+{\cal O}(\omega^{-6}),\end{aligned}$$ where we have used the isotropy of the system and the traceless property of $T^{\alpha\beta}_f$.

Because we consider the $\omega\rightarrow \infty$ limit, we need to take into account the rescaling/mixing effect of the operators. Due to their vanishing anomalous dimensions, the chiral and gluon condensate terms remain unchanged, but the quark component of the energy-momentum tensor changes. To describe its behavior, we rewrite this operator as $$\begin{aligned}
\label{eq:Tf-decompose}
T^{00}_{f}=T'{}^{00}_{f}+\frac{1}{4{C_F}+{N_f}}\left(T^{00}+\frac{2}{{N_f}}\tilde{T}^{00}\right),\end{aligned}$$ where $$\begin{aligned}
T'{}^{00}_{f}\equiv T^{00}_{f}-\frac{1}{{N_f}}\sum_{f'} T^{00}_{f'},\\
T^{00}\equiv \sum_{f'} T^{00}_{f'}+T^{00}_g,\\
\tilde{T}^{00}\equiv 2C_F  \sum_{f'} T^{00}_{f'}-\frac{{N_f}}{2} T^{00}_g .\end{aligned}$$ Here, $T^{\mu\nu}_g\equiv -G^{\mu\alpha}_{a}G^\nu{}_{\alpha a} +g^{\mu\nu}G^2/4$ is the gluon component of the traceless part of the energy-momentum tensor, $N_f$ is the flavor number, and $C_F\equiv ({N_c}^2-1)/(2{N_c})$. We note that $T^{\mu\nu}$ is the traceless part of the full energy-momentum tensor, not the energy-momentum tensor itself. A standard renormalization group (RG) analysis yields the following scaling properties [@Peskin:1995ev]: $$\begin{aligned}
\label{eq:scaling-T'-Ttilde}
\begin{split} 
T'{}^{00}_{f}(\kappa)&=  \left[\frac{\ln\left(\kappa^2_0/\Lambda^2_{\text{QCD}} \right)}{\ln\left(\kappa^2/\Lambda^2_{\text{QCD}} \right)}\right]^{a'} 
T'{}^{00}_{f}(\kappa_0),\\
 \tilde{T}^{00}(\kappa)
&= \left[\frac{\ln\left(\kappa^2_0/\Lambda^2_{\text{QCD}} \right)}{\ln\left(\kappa^2/\Lambda^2_{\text{QCD}} \right)}\right]^{\tilde{a}}
 \tilde{T}^{00}(\kappa_0),
\end{split} \end{aligned}$$ while $T^{00}$ is independent of $\kappa$. Here $\kappa$ and $\kappa_0$ are renormalization scales, $\Lambda_{\text{QCD}}$ is the QCD scale parameter, $a'\equiv 8{C_F}/(3b_0)$, and $\tilde{a}\equiv 2(4{C_F}+{N_f})/(3b_0)$, where $b_0\equiv (11{N_c}-2{N_f})/3$, which appears in the expression $$\begin{aligned}
\label{eq:RG-alphas}
\alpha_s(\kappa)=\frac{4\pi}{b_0 \ln(\kappa^2/\Lambda^2_{\text{QCD}})}.\end{aligned}$$ In the $\omega\rightarrow \infty$ limit, it is natural to choose the RG scale[^3] as $\kappa^2=\omega^2$.

We see that, except for the $T^{00}$ term, all terms in Eq. (\[eq:Tf-decompose\]) are suppressed logarithmically at large $\omega$. Thus, Eqs. (\[eq:OPE-wo-RG-T\]) and (\[eq:OPE-wo-RG-L\]) become $$\begin{aligned}
\nonumber
\delta G_{T}(\omega,{{\mathbf{p}}}) 
&=e^2\sum q^2_f \frac{1}{p^2} \\
\nonumber
&~~~\times \Bigl[\left\{2m_f \delta\left\langle\overline{\psi}_f \psi_f \right\rangle
+\frac{1}{12}\delta\left\langle \frac{\alpha_s}{\pi}G^2\right\rangle \right\}\\
\label{eq:OPE-w-RG-T}
&~~~+\frac{8}{3}\frac{1}{4{C_F}+{N_f}}\frac{\omega^2+{{\mathbf{p}}}^2}{p^2}\delta\left\langle T^{00} \right\rangle 
\Bigr] \nonumber \\
&~~~+{\cal O}(\omega^{-4}),\\
\nonumber
\delta G^{R}_{00}(\omega,{{\mathbf{p}}}) 
&=e^2\sum q^2_f \frac{1}{p^2}\frac{{{\mathbf{p}}}^2}{p^2} \\ 
\nonumber
&~~~\times \Bigl[\left\{2m_f \delta\left\langle\overline{\psi}_f \psi_f \right\rangle
+\frac{1}{12}\delta\left\langle \frac{\alpha_s}{\pi}G^2\right\rangle \right\} \\
\label{eq:OPE-w-RG-L}
&~~~+\frac{8}{3}\frac{1}{4{C_F}+{N_f}} \delta\left\langle T^{00} \right\rangle  
\Bigr]
+{\cal O}(\omega^{-6}). \end{aligned}$$ We note that, in the $\omega\rightarrow \infty$ limit, which is relevant to the derivation of the sum rule, the asymptotic freedom of QCD guarantees that the above expression is exact. In other words, all higher order $\alpha_s$ corrections vanish in this limit.

IR behavior {#ssc:IR}
-----------

On the other hand, the asymptotic behavior in the IR region is described by hydrodynamics [@Kovtun:2012rj], as long as the spatial momentum is small enough. In the channel of the EM current, the basic equations consist of the conservation law and the constituent equation, $$\begin{aligned}
\label{eq:conservation-law}
\partial_0j^0&= -\nabla\cdot {{\mathbf{j}}},\\
\nonumber
{{\mathbf{j}}}&= -D\nabla j^0+\sigma{{\mathbf{E}}}-\sigma\tau_J\partial_0 {{\mathbf{E}}}+\kappa_B\nabla\times{{\mathbf{B}}}\\ 
\label{eq:constituent-eq}
&~~~+{\cal O}(\partial^2 {{\mathbf{E}}}, \partial^2 {{\mathbf{B}}}, \partial^2 j^0 ),\end{aligned}$$ where $D$ is the diffusion constant, $\sigma$ the electrical conductivity and $\tau_J$ and $\kappa_B$ second order transport coefficients corresponding to $\partial_0{{\mathbf{E}}}$ and $\nabla\times{{\mathbf{B}}}$, respectively. ${{\mathbf{E}}}\equiv -\nabla A^0-\partial_0 {{\mathbf{A}}}$ and ${{\mathbf{B}}}\equiv\nabla\times{{\mathbf{A}}}$ are the electric and magnetic fields, where $A^\mu$ is the vector potential.

After performing the Fourier transformation that is defined as $f(p)\equiv \int d^4x e^{ip\cdot x} f(x)$, Eqs. (\[eq:conservation-law\]) and (\[eq:constituent-eq\]) become $$\begin{aligned}
\label{eq:conservation-law-p}
\omega j^0&= {{\mathbf{p}}}\cdot {{\mathbf{j}}},\\
\nonumber
{{\mathbf{j}}}&= -Di{{\mathbf{p}}}j^0+i\sigma (-{{\mathbf{p}}}A^0+\omega {{\mathbf{A}}})\\
\nonumber
&~~~-\sigma\tau_J \omega  (-{{\mathbf{p}}}A^0+\omega {{\mathbf{A}}})-\kappa_B ([{{\mathbf{p}}}\cdot{{\mathbf{A}}}]{{\mathbf{p}}}-{{\mathbf{p}}}^2{{\mathbf{A}}}) \\
\label{eq:constituent-eq-p} 
&~~~+{\cal O}(p^2A, p^2j^0).\end{aligned}$$ Let us solve these equations for the transverse and the longitudinal components of the current. By introducing the transverse component, $j^i_T(p)\equiv P^{ij}_T(p) j^j(p)$, we get the solutions as follows: $$\begin{aligned}
{{\mathbf{j}}}_T(p)&= \left(i\sigma\omega-\sigma\tau_J\omega^2+\kappa_B{{\mathbf{p}}}^2
\right) {{\mathbf{A}}}_T \nonumber \\
&~~~ +{\cal O}(\omega^3{{\mathbf{A}}}_T, \omega{{\mathbf{p}}}^2{{\mathbf{A}}}_T, {{\mathbf{p}}}^4{{\mathbf{A}}}_T),\\
\nonumber
\omega j^0(p)&= -iD{{\mathbf{p}}}^2 j^0
-i{{\mathbf{p}}}^2\sigma A^0 \\
&~~~+ {\cal O}(\omega{{\mathbf{p}}}^2 A^0, {{\mathbf{p}}}^4 A^0, \omega{{\mathbf{p}}}^2 j^0, {{\mathbf{p}}}^4 j^0) \nonumber \\
&~~~+({\text {terms that are proportional to }}{{\mathbf{p}}}\cdot{{\mathbf{A}}}).\end{aligned}$$

By using the linear response theory, the induced current is written as $$\begin{aligned}
\label{eq:linear-response}
j_\mu(p)=-G^R_{\mu\nu}(p)A^\nu(p),\end{aligned}$$ from which we obtain $$\begin{aligned}
\label{eq:hydro-GR-T}
G_{T}(p) &= i\sigma\omega-\sigma\tau_J\omega^2+\kappa_B{{\mathbf{p}}}^2
+{\cal O}(\omega^3, \omega{{\mathbf{p}}}^2, {{\mathbf{p}}}^4),\\
\label{eq:hydro-GR-L}
G^R_{00}(p) &= i\sigma{{\mathbf{p}}}^2\frac{1+{\cal O}(\omega,{{\mathbf{p}}}^2)}{\omega+iD{{\mathbf{p}}}^2+{\cal O}(\omega{{\mathbf{p}}}^2,{{\mathbf{p}}}^4)}.\end{aligned}$$ We note that there is a pole at $\omega=-iD{{\mathbf{p}}}^2$ in the longitudinal channel, which we call the diffusion pole. This is a novel structure that appears only at finite ${{\mathbf{p}}}$. This pole appears as a peak in the spectral function, $$\begin{aligned}
\label{eq:hydro-rho-L}
\rho_{00}(p) &= \sigma{{\mathbf{p}}}^2
\frac{\omega}{\omega^2+(D{{\mathbf{p}}}^2)^2} ,\end{aligned}$$ while at $|{{\mathbf{p}}}|=0$, it reduces to the delta function, $$\begin{aligned}
\label{eq:hydro-rho-L-p=0}
\rho_{00}(\omega,{{\mathbf{0}}}) 
&= \pi\frac{\sigma}{D}\omega\delta(\omega) .\end{aligned}$$

For constructing our sum rules, we, in principle, also need to evaluate the zero temperature part, which should be subtracted later. At $T=0$, Lorentz invariance guarantees the tensor structure of the correlator, $$\begin{aligned}
G^R_{\mu\nu}(p)
&= p^2P^0_{\mu\nu}(p) \tilde{G}^R(p^2),\end{aligned}$$ where $\tilde{G}^R(0)=0$ due to the renormalization condition of the electric charge [@Peskin:1995ev]. The transverse component is given as $G_T(p)= p^2\tilde{G}^R(p^2)$. $\tilde{G}^R(p^2)$ is regular at $p^2=0$ due to the renormalization condition, so it is easy to see that there are no contributions to $\sigma$, $\tau_J$, and $\kappa_B$ at $T=0$, from Eq. (\[eq:hydro-GR-T\]). The $T=0$ contribution only enters in the higher order terms, which are neglected in Eq. (\[eq:hydro-GR-T\]).

On the other hand, the longitudinal component reads $G^R_{00}(p)={{\mathbf{p}}}^2 \tilde{G}^R(p^2)$, so $G^R_{00}(p)$ at $T=0$ and $\omega=0$ is of order ${{\mathbf{p}}}^2 \tilde{G}^R(-{{\mathbf{p}}}^2)$. The $T=0$ contribution does not affect the sum rules 2 (\[eq:sumrule-2-L\]) and 3 (\[eq:sumrule-3-L\]) because as is explained in Sec. \[ssc:sumrule-L\], the relevant quantities to the derivation of these sum rules are $\omega^2\delta G^R_{00}(p)|_{\omega\rightarrow 0}$ and $\omega^4\delta G^R_{00}(p)|_{\omega\rightarrow 0}$, and the $T=0$ contributions vanish. The only sum rule which $T=0$ terms may affect is sum rule 1 (\[eq:sumrule-1-L\]). Nevertheless, we only consider terms in $G^R_{00}(p)|_{\omega\rightarrow 0}$ up to ${{\mathbf{p}}}^2$ in this paper, as is shown in Sec. \[ssc:sumrule-L\]. Therefore, sum rule 1 is also unaffected.

We note that we have so far neglected effects of possible hydro modes, which adds a $\omega^{3/2}$ term at zero momentum to $G^T$ [@Kovtun:2003vj]. It was suggested that this approximation is justified in the large ${N_c}$ limit [@Kovtun:2003vj]. Among our sum rules, sum rule 3 (\[eq:sumrule-3-T\]) in the transverse channel may be changed at finite ${N_c}$ due to this effect. In the analysis of lattice QCD data, this effect practically can be neglected because the IR cutoff of the current lattice QCD is not small enough for such IR energy effects to be detectable.

Sum rules at finite momentum {#sec:sumrule}
============================

In this section, we derive the sum rules at finite momentum in the transverse and the longitudinal channels. For later convenience, we also give the expressions for the sum rules at zero momentum in the transverse channel, which were already obtained in Ref. [@Gubler:2016hnf]. We also confirm that the sum rules are satisfied by the expressions of the spectral function in the chiral and weak coupling limits.

Transverse channel
------------------

Before deriving the sum rules, let us discuss what kind of structure can be expected to appear in the spectral function in a perturbative analysis. First, in the low-energy region, a peak with a width of order $g^4T$ due to the collision effect is expected to appear. This peak is called the transport peak. Its derivation is recapitulated in Appendix \[app:transport-peak\]. At larger energy, $\omega\sim T$, the pair-creation process yields a continuum in the spectral function (see Appendix \[app:continuum\] for its expression). Also, at $\omega\gg T$, the OPE analysis predicts a UV tail, whose derivation is recapitulated in Appendix \[app:UVtail\]. These structures are summarized in Fig. \[fig:spectral-longpaper\], where the expressions in the chiral and weak coupling limits for ${N_c}={N_f}=3$, Eqs. (\[eq:transport-T\]), (\[eq:continuum-T\]), and (\[eq:UVtail-T\]) have been used. The following parameters are chosen for illustrative purposes: $\tau^{-1}/T=0.5$, $|{{\mathbf{p}}}|/T=0.5$, $\kappa_0/T=1$, and $\Lambda_{\text{QCD}}/T=0.67$. We note that the corrections due to the ${{\mathbf{p}}}^2$ terms are almost negligible for this case, though the $|{{\mathbf{p}}}|$ value adopted here is not very small compared to $T$. We furthermore caution that the plots for each structure are reliable only at their energy regions of applicability. Namely, the transport peak is reliable at low energy, the continuum at intermediate and high energy and the UV tail at high energy, respectively. These regions are marked by the vertical lines with attached arrows in the figure. Note that these boundaries are not exact and should only be considered as indicative. One should not take the curves seriously when they are outside of the adequate energy regions.

![ The transport peak, the continuum, and the UV tail of the EM current spectral function in the transverse channel $\delta\rho_T$ as a function of $\omega$. The energy unit is $T$. To draw the figure, the parameters are set as ${N_c}={N_f}=3$, $\tau^{-1}/T=0.5$, $\kappa_0/T=1$, and $\Lambda_{\text{QCD}}/T=0.67$. The solid (dashed) lines correspond to a spatial momentum of $|{{\mathbf{p}}}|/T = 0.5$ ($|{{\mathbf{p}}}| = 0$). The vertical lines with the attached arrows indicate the regions for which the respective analytic expressions can be trusted.[]{data-label="fig:spectral-longpaper"}](spectral-longpaper-T.ranges.eps){width="49.00000%"}

### Sum rule 1

By using the asymptotic expression of $\delta G^R_T$ in the UV and IR energy regions, Eqs. (\[eq:OPE-w-RG-T\]) and (\[eq:hydro-GR-T\]), Eq. (\[eq:derivation\]) for $\mu=\nu=1$ becomes $$\begin{aligned}
\label{eq:sumrule-1-T}
\begin{split}
\kappa_B {{\mathbf{p}}}^2
+{\cal O}({{\mathbf{p}}}^4)
= \frac{2}{\pi}\int^\infty_{0} d\omega
\frac{\delta \rho_{T}(\omega,{{\mathbf{p}}})}{\omega},
\end{split}\end{aligned}$$ where the contribution from the UV part has vanished. This is the first sum rule in the transverse channel (sum rule 1). We note that $|{{\mathbf{p}}}|$ should be small enough to trust this relation, because we have assumed that the IR region is well described by hydrodynamics, which is valid only at small momentum and energy. $\kappa_B$ has been evaluated in lattice QCD [@Brandt:2013faa] with methods that do not suffer from the problem of analytic continuation, so this sum rule can be used to constrain the spectral function. At $|{{\mathbf{p}}}|=0$, it reduces to $$\begin{aligned}
\label{eq:sumrule-1-p0}
\begin{split}
0&= \frac{2}{\pi}\int^\infty_{0} d\omega
\frac{\delta \rho_{T}(\omega,{{\mathbf{0}}})}{\omega},
\end{split}\end{aligned}$$ which was already obtained in Ref. [@Gubler:2016hnf], and also in Ref. [@Bernecker:2011gh] by using the conservation of the EM current.

To get a feeling of how this sum rule is satisfied, let us check the respective contributions from the transport peak, the continuum, and the UV tail below. This also gives us an indication about the sensitivity of the sum rule integral to these three structures. The transport peak in the spectral function at small momentum is given by Eq. (\[eq:transport-T\]), and its contribution to sum rule 1 reads $$\begin{aligned}
\label{eq:T-sumrule1-transport}
\begin{split}
 \frac{2}{\pi}\int^\infty_{0} d\omega
\frac{\delta \rho_{T}(\omega,{{\mathbf{p}}})}{\omega}
&\simeq  {C_{\text{em}}}{N_c}\chi \frac{2}{3} \frac{2}{\pi}
\int^\infty_{0} d\omega \frac{\tau^{-1} }{\omega^2+\tau^{-2}} \\
&~~~\times \left[1 +\frac{{{\mathbf{p}}}^2}{5} \frac{(3\omega^2-\tau^{-2})}{(\omega^2+\tau^{-2})^2}\right]\\
&= {C_{\text{em}}}{N_c}\chi \frac{2}{3},
\end{split}\end{aligned}$$ where ${C_{\text{em}}}\equiv e^2\sum_f q^2_f$, $\tau\sim (g^4T)^{-1}$ is the relaxation time introduced in the Boltzmann equation, and $\chi\equiv T^2/6$. For the leading term we have used $\int^\infty_{0} d\omega \tau^{-1}/(\omega^2+\tau^{-2})=\pi/2$, while the momentum dependent term vanishes because of $\int^\infty_{0} d\omega (3\omega^2-\tau^{-2})/(\omega^2+\tau^{-2})^3=-[\omega/(\omega^2+\tau^{-2})^2]^\infty_0=0$.

The contribution from the continuum (\[eq:continuum-T\]) to sum rule 1 is given as $$\begin{aligned}
\label{eq:T-sumrule1-cont}
\begin{split}
& \frac{2}{\pi}\int^\infty_{0} d\omega
\frac{\delta \rho_{T}(\omega,{{\mathbf{p}}})}{\omega} \\
&\simeq -C_{\text{em}}{N_c}\frac{1}{3\pi^2}
\int^\infty_{0} d\omega \omega
 \Biggl[{n_F}\left(\frac{\omega}{2}\right) \\
&~~~ +{{\mathbf{p}}}^2\left\{\frac{1}{20}{n_F}''\left(\frac{\omega}{2}\right)
 -\frac{1}{\omega^2}{n_F}\left(\frac{\omega}{2}\right)\right\}\Biggr] \\
 &\simeq -C_{\text{em}}{N_c}\frac{1}{3}
 \Biggl[\frac{T^2}{3} 
 -\frac{{{\mathbf{p}}}^2}{2\pi^2}\ln\frac{T}{\mu}\Biggr] ,
\end{split}\end{aligned}$$ where ${n_F}(k^0)\equiv[e^{k^0/T}+1]^{-1}$ is the Fermi distribution function. We have used $\int^\infty_{0} d\omega \omega {n_F}(\omega/2)=\pi^2T^2/3$, and introduced the IR cutoff $\mu$ for the continuum contribution because of a logarithmic IR divergence. The nonsingular parts of the ${{\mathbf{p}}}^2$ terms have been omitted. The result for the continuum (\[eq:continuum-T\]) is obtained from a one-loop calculation, and becomes unreliable when $\omega\lesssim gT$, where the hard-thermal loop resummation becomes necessary [@Moore:2006qn; @Braaten:1990wp]. Therefore, we see that $\mu\sim gT$.

The contribution from the UV tail (\[eq:UVtail-T\]) is estimated as $$\begin{aligned}
\label{eq:T-sumrule1-UVtail}
\begin{split}
 \frac{2}{\pi}\int^\infty_{0} d\omega
\frac{\delta \rho_{T}(\omega,{{\mathbf{p}}})}{\omega}
&\sim {C_{\text{em}}}g^2T^2, {C_{\text{em}}}g^2{{\mathbf{p}}}^2.
\end{split} \end{aligned}$$ Here we have used the fact that the IR cutoff of the UV tail is of order $\sim T$, because the derivation of the UV tail is based on the OPE, which is valid when $\omega\gg T$.

Now let us check that sum rule 1 (\[eq:sumrule-1-T\]) is satisfied. Because $\kappa_B=0$ at weak coupling (see Appendix \[app:transport-peak\]), the integral in Eq. (\[eq:sumrule-1-T\]) needs to vanish in order to satisfy the sum rule. We first see that for the ${{\mathbf{p}}}$-independent part, the contributions from the transport peak and the continuum, which are of order ${C_{\text{em}}}T^2$, cancel while the contribution from the UV tail is of higher order ($\sim {C_{\text{em}}}g^2T^2$), by looking at Eqs.(\[eq:T-sumrule1-transport\]-\[eq:T-sumrule1-UVtail\]). Therefore, the ${{\mathbf{p}}}$-independent part was shown to satisfy sum rule 1 at leading order already in Ref.[@Gubler:2016hnf]. For the ${{\mathbf{p}}}^2$ term, we see that the transport peak does not contribute, and the continuum contribution is of order $\sim {C_{\text{em}}}{{\mathbf{p}}}^2\ln (1/g)$ while the UV tail contribution is suppressed by a factor of $g^2$. Because the continuum contribution is sensitive to the IR cutoff, we need to improve the evaluation by performing the HTL resummation [@Moore:2006qn; @Braaten:1990wp], in order to confirm that this contribution becomes negligible so that sum rule 1 is satisfied at order ${{\mathbf{p}}}^2$. It is furthermore understood that sum rule 1 is mainly sensitive to the transport peak as well as the continuum, while the contribution of the UV tail is small.

### Sum rule 2 {#sssec:sumrule2-T}

In the derivation of Eq.(\[eq:derivation\]), we used only the fact that the retarded Green function is analytic in the upper $\omega$ plane. Thus, we can derive a similar equation in which $\delta G^R(\omega)$ ($\delta \rho(\omega)$) is replaced with $\omega^2\delta G^R(\omega)$ ($\omega^2\delta \rho(\omega)$), $$\begin{aligned}
\begin{split}
&\omega^2\delta G_{T}(\omega,{{\mathbf{p}}})|_{\omega\rightarrow 0}-\omega^2\delta G_{T}(\omega,{{\mathbf{p}}})|_{\omega\rightarrow \infty} \\
&= \frac{2}{\pi}\int^\infty_{0} d\omega \omega
\delta \rho_{T}(\omega,{{\mathbf{p}}}),
\end{split}\end{aligned}$$ for the transverse component. By using Eqs. (\[eq:OPE-w-RG-T\]) and (\[eq:hydro-GR-T\]), this equation becomes $$\begin{aligned}
\label{eq:sumrule-2-T}
\begin{split}
&-e^2\sum q^2_f 
\Biggl[\left\{2m_f \delta\left\langle\overline{\psi}_f \psi_f \right\rangle
+\frac{1}{12}\delta\left\langle \frac{\alpha_s}{\pi}G^2\right\rangle \right\}\\
&~~~+\frac{8}{3}\frac{1}{4{C_F}+{N_f}}\delta\left\langle T^{00} \right\rangle 
\Biggr] 
= \frac{2}{\pi}\int^\infty_{0} d\omega \omega
\delta \rho_{T}(\omega,{{\mathbf{p}}}). 
\end{split}\end{aligned}$$ We call this equation sum rule 2 in what follows. It should be noted that there is no explicit $|{{\mathbf{p}}}|$ dependence on the left-hand side, but again we are implicitly assuming that $|{{\mathbf{p}}}|$ is small enough so that the hydrodynamics well describes the behavior of $\omega^2\delta G_{T}(\omega,{{\mathbf{p}}})|_{\omega\rightarrow 0}$. Also, we emphasize that the expectation values of the local operators on the left-hand side can be evaluated nonperturbatively by lattice QCD without suffering from the problem of analytic continuation. Therefore, sum rule 2 can be used to constrain the shape of the spectral function. As it was already discussed in the previous paper [@Gubler:2016hnf], lattice QCD results show that the left-hand side of Eq.(\[eq:sumrule-2-T\]) is found to be dominated by the $\left\langle T^{00} \right\rangle$ for almost all temperatures around and above $T_c$.

Let us evaluate the contributions to sum rule 2 from the transport peak, continuum, and UV tail. The contribution from the transport peak is found by using Eq. (\[eq:transport-T\]) as $$\begin{aligned}
\label{eq:T-sumrule2-transport}
\begin{split}
\frac{2}{\pi}\int^\infty_{0} d\omega \omega
\delta \rho_{T}(\omega,{{\mathbf{p}}})
&=  {C_{\text{em}}}{N_c}\chi \frac{2}{3} \frac{2}{\pi}\int^\infty_{0} d\omega 
  \frac{\tau^{-1} \omega^2}{\omega^2+\tau^{-2}} \\
&~~~\times \left[1 +\frac{{{\mathbf{p}}}^2}{5} \frac{(3\omega^2-\tau^{-2})}{(\omega^2+\tau^{-2})^2}\right]\\
 &= {C_{\text{em}}}{N_c}\chi \frac{2}{3}
 \left[ \frac{2}{\pi}\tau^{-1}\varLambda +\frac{{{\mathbf{p}}}^2}{5} \right] ,
 \end{split}\end{aligned}$$ where we have introduced the UV cutoff $\varLambda$ for the transport peak because of the linear UV divergence, and used $\int^\infty_{0} d\omega \omega^2(3\omega^2-\tau^{-2})/(\omega^2+\tau^{-2})^3=\pi \tau/2$. As the Boltzmann equation cannot be used when $\omega\gtrsim gT$ since the instantaneous scattering description becomes invalid [@Moore:2006qn], we set $\varLambda\sim gT$.

The contribution from the continuum (\[eq:continuum-T\]) reads $$\begin{aligned}
\label{eq:T-sumrule2-cont}
\begin{split}
& \frac{2}{\pi}\int^\infty_{0} d\omega \omega
\delta \rho_{T}(\omega,{{\mathbf{p}}}) \\
&=  -  C_{\text{em}}{N_c}\frac{1}{3\pi^2} 
\int^\infty_{0} d\omega \omega^3
 \Biggl[{n_F}\left(\frac{\omega}{2}\right) \\
&~~~~ +{{\mathbf{p}}}^2\left\{\frac{1}{20}{n_F}''\left(\frac{\omega}{2}\right)
 -\frac{1}{\omega^2}{n_F}\left(\frac{\omega}{2}\right)\right\}\Biggr] \\
 &= -  C_{\text{em}}{N_c}\frac{1}{45} T^2
 \Biggl[ 14\pi^2T^2  +{{\mathbf{p}}}^2 \Biggr] ,
\end{split}\end{aligned}$$ where we have used $\int^\infty_{0} d\omega \omega^3 {n_F}(\omega/2)=14\pi^4T^4/15$ and $\int^\infty_{0} d\omega \omega^3 {n_F}'' (\omega/2)=8\pi^2T^2$.

The UV tail contribution is, by using Eq. (\[eq:UVtail-T\]), found to be $$\begin{aligned}
\label{eq:T-sumrule2-UVtail}
\begin{split}
 \frac{2}{\pi}\int^\infty_{0} d\omega \omega
\delta \rho_{T}(\omega,{{\mathbf{p}}})
&=  \frac{2}{\pi} {C_{\text{em}}}{N_c}{C_F}\frac{4\pi^2T^4}{27}   \alpha_s(\kappa_0)  \\
&~~~\times \int^\infty_{X_0} dX
 \left[\frac{X_0}{X}\right]^{\tilde{a}+1} \\
&~~~+{\cal O}\left({C_{\text{em}}}g^2T^2{{\mathbf{p}}}^2\right)\\
&=  {C_{\text{em}}}{N_c}{C_F}\frac{4\pi^2T^4}{9}   \frac{2}{4{C_F}+{N_f}}   \\
&~~~+{\cal O}\left({C_{\text{em}}}g^2T^2{{\mathbf{p}}}^2\right),
\end{split}\end{aligned}$$ where we have introduced $X\equiv \ln(\omega/\Lambda_{\text{QCD}})$ and $X_0\equiv \ln(\kappa_0/\Lambda_{\text{QCD}})$, used Eq. (\[eq:RG-alphas\]) and introduced the IR cutoff of the UV tail as $\kappa_0\sim T$.

Let us check whether sum rule 2 is satisfied. From Eqs.(\[eq:T-sumrule2-transport\]-\[eq:T-sumrule2-UVtail\]), we see that for the ${{\mathbf{p}}}$-independent part, the contributions from the continuum and the UV tail have the same order of magnitude ($\sim {C_{\text{em}}}T^4$) while the contribution from the transport peak is much smaller, $\sim {C_{\text{em}}}T^2 \tau^{-1}\varLambda\sim {C_{\text{em}}}g^5 T^4$. These contributions are found to agree with the left-hand side of Eq. (\[eq:sumrule-2-T\]) by using $\delta \langle T^{00}\rangle={N_c}\pi^2T^4(8{C_F}+7{N_f})/60$, which is obtained from Eqs. (\[eq:T00f-free\]) and (\[eq:T00g-free\]). For the ${{\mathbf{p}}}^2$ term, the contributions from the transport peak and the continuum cancel, while the UV tail contribution is much smaller. Therefore, we have confirmed that sum rule 2 is satisfied up to order ${{\mathbf{p}}}^2$, in the chiral and weak coupling limits.

### Sum rule 3

The derivation of the third sum rule turns out to be somewhat more tricky. Equation (\[eq:residue-theorem\]) for $\mu=\nu=1$ can be rewritten as $$\begin{aligned}
\begin{split}
&\delta G_{T}(i\omega,{{\mathbf{p}}})-\delta G_{T}(\infty,{{\mathbf{p}}}) 
= \frac{1}{\pi }\int^\infty_{0} d\omega' \\
&~~~\times\frac{\omega' \delta \rho_{T}(\omega',{{\mathbf{p}}})+\omega[{\mathrm{Re}}\delta G_{T}(\omega',{{\mathbf{p}}})-\delta G_{T}(\infty,{{\mathbf{p}}})]}{\omega'{}^2+\omega^2}\\
&~~~= \frac{2}{\pi }\int^\infty_{0} d\omega' 
\frac{\omega' \delta \rho_{T}(\omega',{{\mathbf{p}}})}{\omega'{}^2+\omega^2} ,
\end{split}\end{aligned}$$ where we have used the relation $$\begin{aligned}
&0 = \nonumber \\
& \int^\infty_{-\infty} d\omega' 
\frac{\omega' \delta \rho_{T}(\omega',{{\mathbf{p}}})-\omega[{\mathrm{Re}}\delta G_{T}(\omega',{{\mathbf{p}}})-\delta G_{T}(\infty,{{\mathbf{p}}})]}{\omega'{}^2+\omega^2},\end{aligned}$$ in the last line, which is obtained by using the residual theorem for the integral $\oint_C d\omega' 
[\delta G_{T}(\omega',{{\mathbf{p}}})-\delta G_{T}(\infty,{{\mathbf{p}}})]/(\omega'+i\omega)$. By subtracting Eq. (\[eq:derivation\]) and $i\omega \delta G_T{}'=i\omega^2\delta G_T{}'2\int^\infty_0 d\omega'/[\pi(\omega^2+\omega'{}^2)]$, which is necessary to regularize the IR singularity in the integral, we get $$\begin{aligned}
\begin{split}
&\delta G_{T}(i\omega,{{\mathbf{p}}})-\delta G_{T}(0,{{\mathbf{p}}}) -i\omega \delta G_T{}'(0,{{\mathbf{p}}}) \\
&= \frac{2}{\pi }\omega^2 \int^\infty_{0} d\omega' \frac{1}{\omega^2+\omega'{}^2}
\left[ \delta \rho_{T}(\omega',{{\mathbf{p}}})\frac{-1}{\omega'} 
+\delta \rho'_T{}(0,{{\mathbf{p}}})\right],
\end{split}\end{aligned}$$ where $'$ stands for the derivative in terms of energy ($\omega$, $\omega'$). Taking the $\omega \to 0$ limit, this reduces to $$\begin{aligned}
\label{eq:derivation-sumrule3-T}
\begin{split}
\frac{1}{2}\delta G_{T}{}''(0,{{\mathbf{p}}})
= \frac{2}{\pi } \int^\infty_{0} d\omega \frac{1}{\omega^3}
\left[ \delta \rho_{T}(\omega,{{\mathbf{p}}})
-\omega \delta \rho'_T{}(0,{{\mathbf{p}}})\right]. 
\end{split}\end{aligned}$$ Here, we have changed the integration variable from $\omega'$ to $\omega$ for simplicity. To get an explicit form of this sum rule, one needs to evaluate $\delta G_{T}{}''(0,{{\mathbf{p}}})$ and $\delta \rho'_T{}(0,{{\mathbf{p}}})$. In the expansion of Eq.(\[eq:hydro-GR-T\]), we get only the $|{{\mathbf{p}}}|=0$ terms as $\delta G_{T}{}''(0,{{\mathbf{p}}})=-2\sigma\tau_J+{\cal O}({{\mathbf{p}}}^2)$ and $\delta \rho'_T{}(0,{{\mathbf{p}}})=\sigma+{\cal O}({{\mathbf{p}}}^2)$. Therefore, we can obtain a sum rule for the $|{{\mathbf{p}}}|=0$ case, which reads $$\begin{aligned}
\label{eq:sumrule-3-T}
-\sigma\tau_J
&= \frac{2}{\pi} \int^\infty_{0} \frac{d\omega}{\omega^3} 
\left[\delta \rho_T(\omega,{{\mathbf{0}}})-\sigma\omega\right].\end{aligned}$$ We call this equation sum rule 3 for the transverse channel. We note that the transport coefficients in the left-hand side cannot be computed by lattice QCD without suffering from the problem of analytic continuation. We do not check that sum rule 3 is satisfied in the chiral and the weak coupling limits since this was already done in our previous paper [@Gubler:2016hnf]. Instead we just cite the order of magnitude of the three contributions: the transport peak contribution is of order ${C_{\text{em}}}T^2 \tau^2\sim {C_{\text{em}}}g^{-8}$ and is equal to the left-hand side of sum rule 3, while the continuum is much smaller, ${C_{\text{em}}}g^{-5}$. The UV tail contribution is the smallest and of order ${C_{\text{em}}}g^{-4}$.

Explicitly taking into account higher order terms in Eq. (\[eq:hydro-GR-T\]) and expanding $\delta\rho_T$ in $|{{\mathbf{p}}}|^2$ order by order, it should be possible to obtain a corresponding sum rule at finite momentum. We leave this task for future work.

Longitudinal channel {#ssc:sumrule-L}
--------------------

Before discussing the sum rules, let us remember that the retarded Green function in the longitudinal channel is exactly known at zero momentum [@Ding:2010ga] from the conservation law of the charge: $$\begin{aligned}
\label{eq:rho-L-p=0-exact}
\rho_{00}(\omega,{{\mathbf{0}}}) 
&= \pi \chi_q \omega\delta(\omega) ,\end{aligned}$$ where $\chi_q\equiv \int d^3{{\mathbf{x}}}\langle j^0({{\mathbf{x}}}) j^0({{\mathbf{0}}}) \rangle/T$ is the charge susceptibility. By matching this result with the hydro result of Eq.(\[eq:hydro-rho-L-p=0\]), we see that the hydro result is exact for all $\omega$ at zero momentum, and $\sigma/D=\chi_q$. For this reason, the sum rules in the longitudinal channel provide nontrivial information only when ${{\mathbf{p}}}$ is finite. Therefore, we consider only the finite momentum case in this subsection.

At finite momentum, the diffusion peak appears in the longitudinal spectral function in addition to the three structures that were already present in the transverse channel, as was explained in Sec. \[ssc:IR\]. To get a feeling about the possible shape of the spectral function in the longitudinal channel, we plot the diffusion peak (\[eq:hydro-rho-L\]), the transport peak (\[eq:transport-00\]), the continuum (\[eq:continuum-00\]), and the UV tail (\[eq:UVtail-00\]) in Fig. \[fig:spectral-longpaper-L\]. The parameters are the same as in Fig. \[fig:spectral-longpaper\]. Again, the approximate regions for which the above analytic descriptions are expected to be valid are indicated by the vertical lines and arrows.

Here we comment on the treatment of the diffusion peak in the traditional QCD sum rule literature: In the conventional sum rule approach, the delta function structure that is similar to the hydro result of Eq.(\[eq:hydro-rho-L-p=0\]) was suggested based on the perturbative calculation in Ref.[@Bochkarev:1985ex], and has been assumed in the subsequent works. Though the two approaches give the same form at $|{{\mathbf{p}}}|=0$ as they should follow the exact results (\[eq:rho-L-p=0-exact\]), the perturbative approach is not generally reliable at $\omega=|{{\mathbf{p}}}|=0$ even when $g$ is small, so that the hydro approach should be adopted. Actually, once we consider finite $|{{\mathbf{p}}}|$, they yield different results.

![ The transport peak, the continuum, and the UV tail of the EM current spectral function in the longitudinal channel $\delta\rho_{00}$ as a function of $\omega$. The energy unit is $T$. The parameters and the meaning of the vertical line with attached arrows are the same as in Fig. \[fig:spectral-longpaper\].[]{data-label="fig:spectral-longpaper-L"}](spectral-longpaper-00.ranges.eps){width="49.00000%"}

### Sum rule 1

The first sum rule for $\delta\rho_{00}$ is obtained from Eq. (\[eq:derivation\]) for $\mu=\nu=0$ by using Eqs. (\[eq:OPE-w-RG-L\]) and (\[eq:hydro-GR-L\]): $$\begin{aligned}
\label{eq:sumrule-1-L}
\begin{split}
\frac{\sigma}{D}
+{\cal O}({{\mathbf{p}}}^2) 
= \frac{2}{\pi}\int^\infty_{0} d\omega
\frac{\delta \rho_{00}(\omega,{{\mathbf{p}}})}{\omega}.
\end{split}\end{aligned}$$ We call this sum rule 1 for the longitudinal channel. Since $\sigma/D$ agrees with the susceptibility $\chi_q$, the left-hand side can be evaluated nonperturbatively by lattice QCD without the problem of analytic continuation.

Let us evaluate the contributions from the diffusion peak, the transport peak, the continuum, and the UV tail, in the weak coupling and the chiral limits. By using Eq. (\[eq:hydro-rho-L\]), the contribution from the diffusion peak is evaluated as $$\begin{aligned}
\begin{split}
\frac{2}{\pi}\int^\infty_{0} d\omega
\frac{\delta \rho_{00}(\omega,{{\mathbf{p}}})}{\omega}
&= \frac{2}{\pi} \sigma{{\mathbf{p}}}^2
\int^\infty_{0} d\omega \frac{1}{\omega^2+(D{{\mathbf{p}}}^2)^2}\\
&= \frac{\sigma}{D},
\end{split}\end{aligned}$$ which is of order ${C_{\text{em}}}T^2$. This contribution is equal to the left-hand side of Eq.(\[eq:sumrule-1-L\]). All the other contributions \[the transport peak (\[eq:transport-00\]), continuum (\[eq:continuum-00\]), and the UV tail (\[eq:UVtail-00\])\] are found to be proportional to ${{\mathbf{p}}}^2$, so the contribution from the diffusion peak is dominant when the momentum is small, and sum rule 1 is satisfied in the limit considered here.

### Sum rule 2 {#sum-rule-2}

The second sum rule is obtained by replacing $\delta G^R(\omega)$ \[$\delta \rho(\omega)$\] with $\omega^2\delta G^R(\omega)$ \[$\omega^2\delta \rho(\omega)$\] in the derivation of sum rule 1, as in Sec. \[sssec:sumrule2-T\]. The result is $$\begin{aligned}
\label{eq:sumrule-2-L}
0= \frac{2}{\pi}\int^\infty_{0} d\omega \omega
\delta \rho_{00}(\omega,{{\mathbf{p}}}).\end{aligned}$$ This is sum rule 2 for the longitudinal channel, which constrains the spectral function. We note that this sum rule can be obtained also by using the current conservation [@Bernecker:2011gh]. Therefore, actually this sum rule is exact, and valid at any value of $|{{\mathbf{p}}}|$, not only for small $|{{\mathbf{p}}}|$.

We again evaluate the contributions to sum rule 2 in the weak coupling and the chiral limits. First we check the diffusion peak contribution. By using Eq. (\[eq:hydro-rho-L\]), we get the order estimate $$\begin{aligned}
\begin{split}
 \frac{2}{\pi}\int d\omega \omega
\delta \rho_{00}(\omega,{{\mathbf{p}}}) 
& =  \frac{2}{\pi}  \sigma{{\mathbf{p}}}^2 
\int^\infty_{0} d\omega \frac{\omega^2}{\omega^2+(D{{\mathbf{p}}}^2)^2} \\
&\sim {C_{\text{em}}}T^2\tau^2 {{\mathbf{p}}}^4 , 
\end{split}\end{aligned}$$ where we have introduced a UV cutoff of the diffusion peak, which is of order $D{{\mathbf{p}}}^2\sim \tau{{\mathbf{p}}}^2$.

The transport peak (\[eq:transport-00\]) contributes as $$\begin{aligned}
\begin{split}
& \frac{2}{\pi}\int^\infty_{0} d\omega \omega
\delta \rho_{00}(\omega,{{\mathbf{p}}}) \\
&=  \frac{2}{\pi} \frac{2}{3}{C_{\text{em}}}{N_c}\chi {{\mathbf{p}}}^2
\int^\infty_{0} d\omega 
 \frac{\tau^{-1}}{\omega^2+\tau^{-2}} \\
&~~~\times \left[1
+{{\mathbf{p}}}^2 \frac{2}{5}\left(\tau^{-2}+\frac{11}{3}\omega^2\right)
\frac{1}{(\omega^2+\tau^{-2})^2}
\right] \\
&=   \frac{2}{3}{C_{\text{em}}}{N_c}\chi {{\mathbf{p}}}^2
 \left[1+\frac{{{\mathbf{p}}}^2  \tau^2}{3} \right] ,
\end{split}\end{aligned}$$ where we have used $\int^\infty_{0} d\omega \tau^{-1}(\tau^{-2}+11\omega^2/3)/(\omega^2+\tau^{-2})^3= \tau^2 (\pi/2)(5/6)$.

The continuum contribution is estimated by using Eq. (\[eq:continuum-00\]), $$\begin{aligned}
\begin{split}
 \frac{2}{\pi}\int^\infty_{0} d\omega \omega
\delta \rho_{00}(\omega,{{\mathbf{p}}})
&= -C_{\text{em}}{N_c}\frac{{{\mathbf{p}}}^2}{3\pi^2} \int^\infty_{0} d\omega \omega \\
&~~~\times \left[{n_F}\left(\frac{\omega}{2}\right)+\frac{{{\mathbf{p}}}^2}{40}{n_F}''\left(\frac{\omega}{2}\right)\right] \\
&= -C_{\text{em}}{N_c}\frac{{{\mathbf{p}}}^2}{3\pi^2}  
 \left[\frac{\pi^2T^2}{3} +\frac{{{\mathbf{p}}}^2}{20}\right], 
\end{split}\end{aligned}$$ where we have used $\int^\infty_{0} d\omega \omega{n_F}''\left(\omega/2\right)=2$.

The UV tail contribution is estimated to be $$\begin{aligned}
\begin{split}
 \frac{2}{\pi}\int^\infty_{0} d\omega \omega
\delta \rho_{00}(\omega,{{\mathbf{p}}})
&= {\cal O}\left( C_{\text{em}} g^2T^2{{\mathbf{p}}}^2, C_{\text{em}} g^2 {{\mathbf{p}}}^4\right),
\end{split}\end{aligned}$$ by using Eq. (\[eq:UVtail-00\]).

By looking at all the contributions, we see that the contributions from the transport peak and the continuum are larger than the one from the UV tail at order ${{\mathbf{p}}}^2$. However, these two contributions are found to cancel each other, so that sum rule 2 is satisfied at leading order in ${{\mathbf{p}}}^2$. For the ${{\mathbf{p}}}^4$ terms, the contributions from the diffusion peak and the transport peak are larger than the ones from the continuum and the UV tail. These contributions are expected to cancel, but we cannot confirm this here since the former contribution was not explicitly calculated.

### Sum rule 3

The third sum rule is obtained by replacing $\delta G^R(\omega)$ \[$\delta \rho(\omega)$\] with $\omega^4\delta G^R(\omega)$ \[$\omega^4\delta \rho(\omega)$\] in the derivation of sum rule 1. It gives $$\begin{aligned}
\label{eq:sumrule-3-L}
\begin{split}
&-e^2\sum q^2_f {{\mathbf{p}}}^2 
\Bigl[\left\{2m_f \delta\left\langle\overline{\psi}_f \psi_f \right\rangle
+\frac{1}{12}\delta\left\langle \frac{\alpha_s}{\pi}G^2\right\rangle \right\} \\
&~~~+\frac{8}{3}\frac{1}{4{C_F}+{N_f}} \delta\left\langle T^{00} \right\rangle  
\Bigr]
= \frac{2}{\pi}\int^\infty_{0} d\omega \omega^3
\delta \rho_{00}(\omega,{{\mathbf{p}}}),
\end{split}\end{aligned}$$ which we call sum rule 3 in the longitudinal channel. We note that there is no ${{\mathbf{p}}}^4$ correction on the left-hand side, but implicitly the smallness of $|{{\mathbf{p}}}|$ is assumed so that hydrodynamics is reliable. As before, the operators on the left-hand side of this sum rule can be evaluated by lattice QCD without having the problem of analytic continuation.

Let us evaluate the contribution to this sum rule in the weak coupling and the chiral limits. The diffusion peak contribution is, by using Eq. (\[eq:hydro-rho-L\]), estimated as $$\begin{aligned}
\begin{split}
 \frac{2}{\pi}\int^\infty_{0} d\omega \omega^3
\delta \rho_{00}(\omega,{{\mathbf{p}}})
&= \frac{2}{\pi}\sigma{{\mathbf{p}}}^2
\int^\infty_{0} d\omega 
\frac{\omega^4}{\omega^2+(D{{\mathbf{p}}}^2)^2}\\
&\sim {C_{\text{em}}}T^2 \tau^4{{\mathbf{p}}}^8,
\end{split}\end{aligned}$$ where we have used the fact that the UV cutoff of the diffusion peak is of order $\tau {{\mathbf{p}}}^2$.

The transport peak contributes as $$\begin{aligned}
\begin{split}
& \frac{2}{\pi}\int^\infty_{0} d\omega \omega^3
\delta \rho_{00}(\omega,{{\mathbf{p}}}) \\
&=  \frac{2}{\pi}\frac{2}{3}{C_{\text{em}}}{N_c}\chi {{\mathbf{p}}}^2 
\int^\infty_{0} d\omega  \frac{\omega^2 \tau^{-1}}{\omega^2+\tau^{-2}} \\
&~~~\times \left[1
+ \frac{2}{5}\left(\tau^{-2}+\frac{11}{3}\omega^2\right)
\frac{{{\mathbf{p}}}^2}{(\omega^2+\tau^{-2})^2} \right] \\
&= {\cal O}({C_{\text{em}}}T^2{{\mathbf{p}}}^2\tau^{-1}\varLambda)
+ {C_{\text{em}}}{N_c}\chi {{\mathbf{p}}}^4 \frac{2}{5} ,
\end{split}\end{aligned}$$ where we have used $\int^\infty_{0} d\omega \left(\tau^{-2}+11\omega^2/3\right) \omega^2 \tau^{-1}/(\omega^2+\tau^{-2})^3=(\pi/2)(3/2)$.

The contribution from the continuum (\[eq:continuum-00\]) reads $$\begin{aligned}
\begin{split}
 \frac{2}{\pi}\int^\infty_{0} d\omega \omega^3
\delta \rho_{00}(\omega,{{\mathbf{p}}})
&= -C_{\text{em}}{N_c}\frac{{{\mathbf{p}}}^2}{6\pi}  \frac{2}{\pi}\int^\infty_{0} d\omega \omega^3 \\
&~~~\times \left[{n_F}\left(\frac{\omega}{2}\right)+\frac{{{\mathbf{p}}}^2}{40}{n_F}''\left(\frac{\omega}{2}\right)\right] \\
&= -C_{\text{em}}{N_c}\frac{{{\mathbf{p}}}^2}{3}  \\
&~~~\times \left[\frac{14\pi^2T^4}{15}+ \frac{{{\mathbf{p}}}^2}{5}T^2\right]. 
\end{split}\end{aligned}$$

Finally, the UV tail contribution (\[eq:UVtail-00\]) is $$\begin{aligned}
\begin{split}
 \frac{2}{\pi}\int^\infty_{0} d\omega \omega^3
\delta \rho_{00}(\omega,{{\mathbf{p}}})
&= \frac{2}{\pi}{C_{\text{em}}}\alpha_s(\kappa_0) {N_c}{C_F}\frac{4\pi^2T^4}{27} {{\mathbf{p}}}^2 \\
&~~~\times \int^\infty_{X_0} dX
  \left[\frac{X_0}{X}\right]^{\tilde{a}+1} \\ 
&~~~+{\cal O}({C_{\text{em}}}g^2T^2{{\mathbf{p}}}^4)\\
  &= {C_{\text{em}}}{N_c}{{\mathbf{p}}}^2 \frac{8\pi^2T^4}{9} 
 \frac{{C_F}}{4{C_F}+{N_f}} \\
&~~~ +{\cal O}({C_{\text{em}}}g^2T^2{{\mathbf{p}}}^4).
\end{split}\end{aligned}$$

Let us check the sum rule order by order in ${{\mathbf{p}}}^2$. For the ${{\mathbf{p}}}^2$ terms, we see that the contributions from the continuum and the UV tail are dominating, and their sum is equal to $-C_{\text{em}}{N_c}{{\mathbf{p}}}^2\pi^2T^4 2(8{C_F}+7{N_f})/[45(4{C_F}+{N_f})]$, which agrees with the left-hand side of sum rule 3 (\[eq:sumrule-3-L\]). The order ${{\mathbf{p}}}^4$ terms are dominated by the transport peak and the continuum, but they cancel each other. Therefore, sum rule 3 is shown to be satisfied up to order ${{\mathbf{p}}}^4$.

### Sum rules for $\delta\rho_L$

Before ending this section, we note that two sum rules for $\delta\rho_L$ can be derived from our three sum rules for $\delta\rho_{00}$, (\[eq:sumrule-1-L\]), (\[eq:sumrule-2-L\]), and (\[eq:sumrule-3-L\]). Using $\rho_L(p)=p^2 \rho_{00}(p)/{{\mathbf{p}}}^2$, we get $$\begin{aligned}
& -\frac{\sigma}{D}
+{\cal O}({{\mathbf{p}}}^2) 
= \frac{2}{\pi}\int^\infty_{0} d\omega
\frac{\delta \rho_{L}(\omega,{{\mathbf{p}}})}{\omega}, \\
\nonumber
& -e^2\sum q^2_f  
\Bigl[\left\{ 2m_f  \delta\left\langle\overline{\psi}_f \psi_f \right\rangle
+\frac{1}{12}\delta\left\langle \frac{\alpha_s}{\pi}G^2\right\rangle \right\} \\
&+\frac{8}{3}\frac{1}{4{C_F}+{N_f}} \delta\left\langle T^{00} \right\rangle  
\Bigr] 
= \frac{2}{\pi}\int^\infty_{0} d\omega \omega
\delta \rho_{L}(\omega,{{\mathbf{p}}})\end{aligned}$$ Naively, one would expect these to agree at $|{{\mathbf{p}}}|=0$ with the corresponding sum rules for $\delta \rho_T$, Eqs.(\[eq:sumrule-1-T\]) and (\[eq:sumrule-2-T\]). This is indeed the case for the latter sum rule, but the former one does not agree, because of the $-\sigma/D$ term on its left-hand side. This disagreement can be traced back to the singularity of $\delta \rho_L$ at the diffusion pole $\omega=-iD{{\mathbf{p}}}^2$.

Application to lattice QCD data analysis {#sec:lattice}
========================================

In this section, we demonstrate that our sum rules can be used to improve a fit of the spectral function to the latest Euclidean time lattice QCD data given in Ref.[@Brandt:2015aqk]. In doing this, we make use of the sum rules in three different ways.

1.  Providing guidance in the choice of the functional forms used to parametrize the spectral function.

2.  Reducing the number of fitting parameters (sum rules 1 and 2).

3.  Determining transport coefficients (sum rule 3).

We note that we use all the sum rules including sum rule 2, which was missing in the previous work [@Gubler:2016hnf].

In the following discussion of lattice QCD data, we restrict ourselves to the zero-momentum case, as most currently available lattice QCD data are provided only in this limit. This simplifies the situation in the sense that for $|{{\mathbf{p}}}| = 0$ the correlator has only one independent component and some parameters drop out of the sum rules \[such as $\kappa_B$ in Eq.(\[eq:sumrule-1-T\])\]. In this section, we use the notation $\rho(\omega, T) = \rho_{T}(\omega, |{{\mathbf{p}}}|=0) = \rho_{L}(\omega, |{{\mathbf{p}}}|=0)$ for the spectral function. The application of our sum rules to nonzero momentum lattice QCD data is left for future work. We hence only use the sum rules of Eqs.(\[eq:sumrule-1-p0\]) and (\[eq:sumrule-2-T\]) at $|{{\mathbf{p}}}| = 0$ and (\[eq:sumrule-3-T\]).

Parametrization of the spectral function in vacuum and at finite temperature
----------------------------------------------------------------------------

Our goal here is to find a functional form of the spectral function that is consistent with the sum rules 1-3, that is, that does not lead to divergent results for these sum rules. Note that some terms of the finite temperature spectra used in Ref.[@Brandt:2015aqk] in fact lead to divergences for sum rules 2 and 3 and therefore violate them.

First, we start with the spectral function in vacuum, for which we can follow the parametrization of Ref.[@Brandt:2015aqk], $$\begin{aligned}
& \frac{1}{C_{\mathrm{em}}} \rho(\omega, T \simeq  0) \nonumber \\
=& \,\, \frac{\pi}{3} a_V \delta(\omega - m_V)  \nonumber \\
&\,\, + (1 + k_1) \frac{1}{4 \pi} \Theta(\omega - \Omega_0) 
                                                \omega^2 \tanh \Bigl( \frac{\omega \beta_0}{4} \Bigr)
\label{eq:vac}                                                \end{aligned}$$ This form is adapted to our notation, which differs from that of Ref.[@Brandt:2015aqk] by a factor of $1/6$ and the treatment of $C_{\mathrm{em}}$. Here $\beta_0 = 1/T_0$, where $T_0$ corresponds to the low temperature of the “vacuum" lattice ensemble of Ref.[@Brandt:2015aqk] ($T_0 \simeq 32\,\mathrm{MeV}$). The values of the parameters $a_V$, $m_V$, $k$ and $\Omega_0$ are determined by the fit. Here, the $\delta$-function peak corresponds to the $\rho$ meson, as an isospin 1 current was used in Ref.[@Brandt:2015aqk]. To check whether the above is a reasonable parametrization, we have performed a MEM analysis [@Jarrell:1996rrw; @Asakawa:2000tr; @Gubler:2010cf] of the vacuum data provided in Ref.[@Brandt:2015aqk]. In this simple analysis, we have ignored correlations between data points at different time slices. The result is shown as a red line in Fig.\[fig:MEM.fit\]. Here, we have chosen the default model, which is an input of the MEM algorithm, to match the perturbative value of the spectral function at high energy \[$\rho(\omega)/\omega^2 =  3/(2\pi)$, blue dashed line in Fig.\[fig:MEM.fit\]\].

![Result of a MEM analysis of vacuum data of Ref.[@Brandt:2015aqk] (red solid line), compared with a simple fit using Eq.(\[eq:vac\]) (black solid line). The blue dashed line shows the default model used in the MEM analysis.[]{data-label="fig:MEM.fit"}](MEM.fit.3.eps){width="49.00000%"}

It is seen in the figure that at low energy the spectral function is dominated by a single peak, while at high energy, there is an almost flat continuum. We furthermore see hints of excited states at around $a\omega \simeq 0.8$, but their effect seems to to be weak. In all, Eq.(\[eq:vac\]) turns out to be a reasonable, if somewhat rough, parametrization of the spectral function at temperatures much below $T_c$. As is also shown in Fig.\[fig:MEM.fit\] and further discussed later, the fit results for the peak position $m_V$ and the onset of the continuum $\Omega_0$ agree well with those obtained from the MEM analysis.

Next, we consider the spectral function at temperatures around and above $T_c$, where large modifications are expected. As it is already discussed in detail in this work, we most importantly expect the appearance of a transport peak at $\omega \simeq 0$ and an UV tail at high energy. We again follow the parametrization used in Ref.[@Brandt:2015aqk], but modify it such that sum rules 1-3 can be satisfied. Specifically, we use $$\begin{aligned}
& \frac{1}{C_{\mathrm{em}}} \rho(\omega, T) \nonumber \\
=&\,\, \frac{\omega A_T \Gamma_T}{3(\Gamma_T^2 + \omega^2)} \bigl[ 1 - A(\omega) \bigr]
                                + \frac{\pi}{3} a_T \delta(\omega - m_V) \nonumber \\ 
                                 &\,\, + A(\omega) \bigl[1 + \tilde{k}(\omega) \bigr] \frac{1}{4 \pi} 
                                 \omega^2 \tanh \Bigl( \frac{\omega \beta_T}{4} \Bigr) \nonumber \\
                                &\,\, + \frac{c_0}{4 \pi} \Theta(\omega - \Omega_0) \frac{1}{\omega^2} \frac{1}{\bigl[\ln(\omega/\Lambda_{\mathrm{QCD}})\bigr]^{1+ \tilde{a}}},  
\label{eq:finiteT}                                                \end{aligned}$$ with $$\begin{aligned}
A(\omega) = \tanh\Bigl( \frac{\omega^2}{\Delta^2}\Bigr) 
\label{eq:A}\end{aligned}$$ and $$\begin{aligned}
\tilde{k}(\omega) = k_1 + k_2 \Biggl[1 - \tanh^2\Bigl( \frac{\omega}{\Omega_0 \eta} \Bigr) \Biggr]. 
\label{eq:k}\end{aligned}$$ Let us here mention the differences between our above parametrization and that of Ref.[@Brandt:2015aqk]. First, the factors $\bigl[ 1 - A(\omega) \bigr]$ and $A(\omega)$, which were introduced already in our previous work [@Gubler:2016hnf], are there for cutting off the divergence in the integral of sum rules 2 and 3. Furthermore, the factor $1/\bigl[\ln(\omega/\Lambda_{\mathrm{QCD}})\bigr]^{1+ \tilde{a}}$, which is adapted from the perturbative expression, is employed to make integral of sum rule 2 finite.

Lattice data used in the fit
----------------------------

Reference [@Brandt:2015aqk] provides correlator data of altogether five lattice ensembles with respective temperatures $32$, $169$, $203$, $254$ and $338\,\mathrm{MeV}$, which in the setting of that work corresponds to $0.15$, $0.8$, $1.0$, $1.25$ and $1.76\,T_c$. The spatial extent of the lattice is $N_s = 64$ for all ensembles, while the temporal size is $N_{\tau} = 128$, 24, 20, 16 and 12. The lattice spacing is $a = 0.0486\,\mathrm{fm}$ and the vacuum pion mass $270\,\mathrm{MeV}$, which shows that this is not yet a physical point simulation. For more details, we refer the reader to the original publication of Ref.[@Brandt:2015aqk].

The lattice data are related to the spectral function $\rho(\omega)$ through the following integral: $$\begin{aligned}
G^{\mathrm{E}}(\tau, T) = \int_{0}^{\infty} \frac{d \omega}{2 \pi} \frac{\rho(\omega,T)}{C_{\mathrm{em}}} \frac{\cosh[\omega(\tau - 1/2T)]}{\sinh(\omega/2T)}.
\label{eq:Eucl.vec.correlator}\end{aligned}$$ Here, $\tau$ is defined in the interval $0 \leq \tau < 1/T$ and is symmetric with respect to the central point $\tau = 1/(2T)$. Hence, only half of the temporal data points can be used as independent information for the fit. Furthermore, for small $\tau$ values, the data could contain lattice artifacts. We follow Ref.[@Brandt:2015aqk] and use data points $\tau/a \in [4:48]$ for the vacuum ensemble and $\tau/a \in [4:N_{\tau}/2]$ for the others.

Treatment of the parameters
---------------------------

We treat the various parameters introduced in Eqs.(\[eq:vac\])-(\[eq:k\]) as follows. We keep $m_V$, $k_1$ and $\Omega_0$ fixed for both vacuum and finite temperature. $\eta$, which only enters in the nonzero temperature parametrization, is also kept fixed for all ensembles. All other parameters, $A_T$, $\Gamma_T$, $\Delta$, $a_T$, $k_2$ and $c_0$ are allowed to depend on temperature. To reduce the numbers of unknowns to be fitted, we however make use of the sum rules 1 and 2 to constrain the parameters. To be specific, the sum rules are used to fix $A_T/\Gamma_T$ and $c_0$ at each temperature. To be able to use sum rule 2, one, in principle, needs the quark- and gluon-condensate values and the energy density as a function of temperature. For the simple trial analysis of this work, we employ the values provided in Ref.[@Bazavov:2014pvz], which we rescale to two flavors according to the ratio in the perturbative high-temperature limit. As the quark masses and other lattice parameters in Ref.[@Bazavov:2014pvz] are different from those of Ref.[@Brandt:2015aqk], this is not a completely consistent treatment and will have to be improved in a more quantitative analysis in the future.

Fit results
-----------

Performing our fit, it turns out that even if one uses both sum rules 1 and 2 as constraints, the lattice data are not sufficient to determine all parameters of Eqs.(\[eq:finiteT\])-(\[eq:k\]) with good accuracy. We found for instance that certain combinations of the ratio $A_T/\Gamma_T$, which is proportional to the electric conductivity $\sigma_{\mathrm{el}}$, and $\Delta$ lead to similar $\chi^2$ values. This means that the systematic uncertainty of these quantities is large and that more data points with smaller errors will be needed to further constrain the parameters, especially for the lattice ensembles with temperatures above $T_c$, where the number of usable data points is rather small. In the following we provide one possible fit result, which is obtained by demanding that $\sigma_{\mathrm{el}}/T$ is a monotonically increasing function of $T$ and take values in reasonable agreement with previous works [@Brandt:2015aqk; @Aarts:2007wj; @Francis:2011bt; @Ding:2010ga; @Ding:2013qw; @Ding:2014dua; @Amato:2013naa; @Brandt:2012jc]. The numerical result of the fit is given in Table\[tab:fit.values\]. These give a $\chi^2$ value of $\chi^2/\mathrm{d.o.f.} = 0.87$. We emphasize once more that this is not the only solution and should therefore be regarded as an illustrative example rather than the final result. Nevertheless, the purpose of our analysis is to demonstrate the possibility to use the sum rules to improve the lattice QCD analysis, and to this end, we believe that our fit suffices.

   Fitted parameters     Fixed     $0.8\,T_c$   $1.0\,T_c$   $1.25\,T_c$   $1.67\,T_c$
  -------------------- ---------- ------------ ------------ ------------- -------------
         $m_V$           0.206         -            -             -             -
         $a_V$          0.000860       -            -             -             -
       $\Omega_0$        0.319         -            -             -             -
         $k_1$           0.0909        -            -             -             -
         $\eta$           2.04         -            -             -             -
       $a_T/T^3$           -          4.15         2.94         2.15          1.08
      $\Gamma_T/T$         -          1.36         0.54       $\infty$      $\infty$
   $A_T/(\Gamma_T T)$      -         0.450        0.728         0.728         0.758
      $\Delta_T/T$         -          4.06         2.97         0.79          1.23
       $c_0/T^3$           -          -166         -250         -186          - 95
         $k_2$             -        -0.0294      -0.0377       -0.0322       -0.0040

Let us now discuss the resulting spectral functions in some detail. First, we compare in Fig.\[fig:MEM.fit\] the vacuum spectral function of Eq.(\[eq:vac\]) with that obtained by MEM. It is seen that the two methods indeed give a quantitatively similar result, while the details naturally disagree due to the roughness of the parametrization of Eq.(\[eq:vac\]) and the limited resolution of MEM. Next, we consider the spectral functions for temperatures around and above $T_c$. The spectra are shown in Fig.\[fig:full.spec\].

![The fitted spectral functions for temperatures between $0.8$ and $1.67\,T_c$. To improve the visibility, we have transformed the $\rho$ meson peak, which in the fit is assumed to be a delta function, into a Gaussian with a small, but nonzero width.[]{data-label="fig:full.spec"}](full.spec.2.eps){width="49.00000%"}

It is seen in this figure that the transport peak gradually grows with increasing temperature. One furthermore observes from the behavior of the discontinuity around $a\omega\simeq 0.3$ that the UV tail, which is parametrized to appear above $\Omega_0$, grows with larger temperature. We note that the coefficient of the UV tail ($c_0$) is negative in this fit, which might indicate that $\langle \tilde{T}^{00}(\omega\sim T)\rangle $ is negative, as is suggested from Eq. (\[eq:deltarho\_T\_UVtail\]). Also note that we here have divided the spectral function by $\omega$, which means that it approaches a constant at small $\omega$, while it approaches a form proportional to $\omega$ at large energy.

To examine the behavior of the spectral function more in detail, we show the same spectral functions with the vacuum part subtracted in Fig.\[fig:subtr.spec\].

![Same as Fig.\[fig:full.spec\], but with vacuum part subtracted.[]{data-label="fig:subtr.spec"}](subtr.spec.2.eps){width="49.00000%"}

At low energy, the same features can be observed as in the previous figure. At high energy, the dominant part of the continuum gets subtracted and only the UV tail remains. Also note that because the residue of the $\rho$ meson mass is reduced at finite temperature, the residue becomes negative for the subtracted spectral function shown here.

Next, let us examine the quality of our fit by numerically integrating the spectral function as in Eq.(\[eq:Eucl.vec.correlator\]) and compare the result with the lattice data of Ref.[@Brandt:2015aqk]. This comparison is shown as red and blue data points in Fig.\[fig:comp.latt\].

![image](T0.0.normalized.comp.forpaper3.eps){width="45.00000%"} ![image](T0.8.normalized.comp.forpaper1.eps){width="45.00000%"} ![image](T1.0.normalized.comp.forpaper.eps){width="45.00000%"} ![image](T1.25.normalized.comp.forpaper.eps){width="45.00000%"} ![image](T1.67.normalized.comp.forpaper.eps){width="45.00000%"}

In these plots we show the ratio between the lattice data and the integrated spectral function, which should be 1 for a perfect fit. Overall, it can be seen that both our and the fit of Ref.[@Brandt:2015aqk] can reproduce the lattice data fairly well, especially if one considers the very small errors of most of the lattice data points. For the lowest temperature ($T = 0.15\,T_c$), both fits exhibit almost the same behavior for small $\tau/a$ values, while for $\tau/a \gtrsim 40$, the fit of Ref.[@Brandt:2015aqk] works slightly better than ours. At higher temperatures, both fits perform equally well, being consistent with the lattice results within errors for most $\tau/a$ values. This shows that the available lattice data are not sufficient to distinguish our finite temperature spectral functions from those obtained in Ref.[@Brandt:2015aqk].

Comparison with fit of Ref.[@Brandt:2015aqk]
--------------------------------------------

In Ref.[@Brandt:2015aqk] a fit similar to ours was performed with the following functional form. For the vacuum part, the same parametrization as in Eq.(\[eq:vac\]) was employed, while for nonzero temperature $$\begin{aligned}
& \frac{1}{C_{\mathrm{em}}} \rho(\omega, T) \nonumber \\
&= \,\, \frac{\omega A_T \Gamma_T}{3(\Gamma_T^2 + \omega^2)} 
                                + \frac{\pi}{3} a_T \delta(\omega - m_V)  \nonumber \\
&~~~ + \Theta(\omega - \Omega_T) \bigl[1 + \tilde{k}(\omega) \bigr] \frac{1}{4 \pi} 
                                 \omega^2 \tanh \Bigl( \frac{\omega \beta_T}{4} \Bigr) \nonumber \\
&~~~ + \frac{c_0}{4 \pi} \Theta(\omega - \Omega_0) \frac{1}{\omega^2} 
\label{eq:finiteT.Brandt}                                                \end{aligned}$$ was used, which we have again adapted to our notation. It should be noted that this form is contradicting with the sum rule 2, as its first and fourth terms lead to divergences in that sum rule. If one however only uses sum rule 1, a fit is possible. The authors of Ref.[@Brandt:2015aqk] have tested several models, in which some combinations of parameters are set to 0 or other fixed values for certain temperatures. For details, we refer the reader to Ref.[@Brandt:2015aqk].

Here, we compare our spectral function with the results obtained for model 2c in Ref.[@Brandt:2015aqk]. The other models of that work have different features in certain regions of $\omega$, but the same overall behavior. In Fig.\[fig:full.spec.comp\] we show the full (nonsubtracted) spectral functions for four different temperatures.

![image](full.spec.comp.T.0.8.forpaper.2.eps){width="45.00000%"} ![image](full.spec.comp.T.1.0.forpaper.2.eps){width="45.00000%"} ![image](full.spec.comp.T.1.25.forpaper.2.eps){width="45.00000%"} ![image](full.spec.comp.T.1.67.forpaper.2.eps){width="45.00000%"}

It is seen in these figures that while the details of model 2c differ from our spectral function, the general structure is the same. The biggest difference between our spectral function and those of Ref.[@Brandt:2015aqk] can be found for the $T = 0.8\,T_c$ case in the region around the $\rho$ meson peak, where all models of Ref.[@Brandt:2015aqk] are close to 0, while our spectral function is smoothly connected to the continuum. The true spectral function in this region likely lies between these two extremes, as the $\rho$ meson in reality has quite a large width and is placed on top of a smooth $\pi \pi$ continuum, but still is the dominating structure below energies of about $1$ GeV (see, for instance, Fig.1 of Ref.[@Kwon:2008vq]).

Evaluation of physical quantities
---------------------------------

From the above fit results, we can determine the electric conductivity from the numerical values of $A_{T}/\Gamma_{T}$. In our parametrization, it is given as $$\begin{aligned}
\frac{\sigma_{\mathrm{el}}}{C_{\mathrm{em}}} = \frac{1}{C_{\mathrm{em}}} \lim_{\omega \to 0} \frac{\rho(\omega,T)}{\omega} = \frac{1}{3} \frac{A_T}{\Gamma_T}.  
\label{eq:elec.conduct}\end{aligned}$$ From the values of Table\[tab:fit.values\], $\sigma_{\mathrm{el}}$ is obtained as shown in the second column of Table\[tab:le.constants\] and on the upper panel of Fig.\[fig:low.en.const\].

![The electric conductivity $\sigma_{\mathrm{el}}$ and the second order transport coefficient $\tau_{J}$ as functions of the temperature $T$. Both quantities are made dimensionless by multiplying appropriate powers of $T$. For $\sigma_{\mathrm{el}}$, the values obtained in Ref.[@Brandt:2015aqk] are also shown. Here, we have for simplicity combined the various errors quoted in Ref.[@Brandt:2015aqk] into a single one for each data point.[]{data-label="fig:low.en.const"}](sigma.2.eps "fig:"){width="49.00000%"} ![The electric conductivity $\sigma_{\mathrm{el}}$ and the second order transport coefficient $\tau_{J}$ as functions of the temperature $T$. Both quantities are made dimensionless by multiplying appropriate powers of $T$. For $\sigma_{\mathrm{el}}$, the values obtained in Ref.[@Brandt:2015aqk] are also shown. Here, we have for simplicity combined the various errors quoted in Ref.[@Brandt:2015aqk] into a single one for each data point.[]{data-label="fig:low.en.const"}](tau.eps "fig:"){width="49.00000%"}

We however emphasize once again here, that other fits with similar $\chi^2$ values but rather different electric conductivities are possible and that the numbers shown here just represent one of many possible solutions. On the upper panel of Fig.\[fig:low.en.const\], we furthermore show the results obtained in Ref.[@Brandt:2015aqk] for comparison. It is seen that with the exception of the point at $T = 1.67\,T_c$, the results are not consistent, even though they show the same general tendency. This once more indicates that the systematic uncertainty in the evaluation of this quantity is still rather large.

Next, we can now use sum rule 3 to estimate the value of the second order transport coefficient $\tau_{J}$, as all other ingredients in that sum rule are known. The results of such a computation are found in the third column of Table\[tab:le.constants\] and on the lower plot of Fig.\[fig:low.en.const\].

       $T$       $\sigma_{\mathrm{el}}/(C_{\mathrm{em}}T)$   $T\tau_{J}$
  ------------- ------------------------------------------- -------------
   $0.8\,T_c$                      0.150                        1.21
   $1.0\,T_c$                      0.243                        2.16
   $1.25\,T_c$                     0.243                        1.96
   $1.67\,T_c$                     0.253                        1.29

  : The electric conductivity $\sigma_{\mathrm{el}}$ and the second order transport coefficient $\tau_{J}$ at various temperatures, as obtained from the fit to lattice QCD data.[]{data-label="tab:le.constants"}

It is seen in this figure that $\tau_{J}$ exhibits quite an interesting behavior as a function of $T$. Namely, it increases for temperatures below $T_c$, takes a maximum at $T = T_c$ and then decreases again for temperatures above $T_c$. It remains to be seen whether this behavior is an artifact of our fit and/or our parametrization of the spectral function or if it is a real physical effect.

As a last result, we give the thermal dilepton rate $dN_{l^{+}l^{-}}/d\omega d^3 {{\mathbf{p}}}$ for vanishing spatial momentum ($|{{\mathbf{p}}}| = 0$), which can be easily obtained from the relation between the spectral function and the dilepton rate: $$\begin{aligned}
\frac{dN_{l^{+}l^{-}}}{d\omega d^3 {{\mathbf{p}}}}(|{{\mathbf{p}}}|=0)  = \frac{\alpha_{\mathrm{em}}^2}{\pi^3 \omega^2} 
\frac{\rho(\omega,T)}{e^{\omega/T} - 1}. \end{aligned}$$ The results are shown in Fig.\[fig:dilepton.rate\], where we have adjusted the horizontal axis to physical units (MeV) and where we also show the corresponding model 2c results of Ref.[@Brandt:2015aqk]. We see that our result is larger than the one in Ref.[@Brandt:2015aqk], especially at $T=0.8T_c$. This difference can be understood from the absence of the spectrum near the vector meson peak in Ref.[@Brandt:2015aqk], which was discussed in the previous subsection.

![The thermal dilepton rate for two-flavor QCD, extracted from our fit to lattice QCD data of Ref.[@Brandt:2015aqk]. The thick solid lines show the result of our fit, while the thin dashed lines correspond to model 2c of Ref.[@Brandt:2015aqk].[]{data-label="fig:dilepton.rate"}](dilepton.mod.2c.eps){width="49.00000%"}

Summary and Concluding Remarks {#sec:summary}
==============================

In the first part of this paper, we derived and discussed five exact sum rules (two for the transverse and three for the longitudinal parts) for the vacuum-subtracted spectral functions of the vector channel at finite temperature, which are determined completely by the UV and IR behavior of the vector correlator. The UV part can be obtained from the OPE, while the IR behavior is accurately described by hydrodynamics. The sum rules are valid for nonzero momentum ${{\mathbf{p}}}$, which should however be small enough such that the hydrodynamic description of the vector correlator in the IR regime can be trusted. In the limit $|{{\mathbf{p}}}| \to 0$, transverse and longitudinal parts approach each other such that the sum rules which we have already derived in Ref.[@Gubler:2016hnf] remain.

In perturbation theory, it has been known that three (four) distinct structures emerge in the transverse (longitudinal) channel: a transport peak, an exponentially suppressed continuum, and a power suppressed UV tail in both channels, and a diffusion peak in the longitudinal channel. All these structures can, in principle, contribute to the various sum rules and therefore need to be taken into account to test their validity. Doing this, we found that the sum rules are indeed satisfied in the weak coupling regime. This exercise also gives a rough idea about how the different parts of the full spectral function can be expected to contribute in different ways to each sum rule.

In the second part of the paper, we have employed recent two-flavor dynamical lattice QCD data at almost zero and finite temperature and at zero momentum to perform a spectral fit, in which two sum rules (1 and 2) are used as constraints to reduce the number of parameters to be fitted. The third sum rule (3) in turn enables us to extract the value of the second order transport coefficient $\tau_J$. We note that the second sum rule was not used in our previous work [@Gubler:2016hnf]. It, however, needs to be emphasized here that we found the fit not to be completely stable in the sense that we confirmed the existence of several local minima with comparable $\chi^2$ values, which means that more data points with increased precision will be needed for uniquely determining the true shape of the spectral function. Nevertheless, we succeeded in demonstrating that the sum rules can be used to improve the lattice QCD data analysis.

In our fit, we employed lattice data at zero momentum and have therefore used the three sum rules already derived in Ref.[@Gubler:2016hnf]. Once lattice data at nonzero momentum are available, it would be interesting to apply the sum rules derived in this paper to their analysis. Also, once the lattice QCD analysis at the physical point and in the continuum limit becomes available, it will be possible to use the phenomenological form of the vector spectral function at $T=0$ obtained from experiment. We leave these topics for future work.

Acknowledgements {#acknowledgements .unnumbered}
================

The research of D. S. is supported by the Alexander von Humboldt Foundation. The research of P.G. is supported by the Mext-Supported Program for the Strategic Foundation at Private Universities, “Topological Science" (No. S1511006). We thank Bastian B. Brandt for fruitful discussions.

Evaluation of the transport peak at weak coupling {#app:transport-peak}
=================================================

In this appendix, we evaluate the transport peak appearing in the spectral function by using the Boltzmann equation in the relaxation time approximation for massless quarks. This appendix is in essence a recapitulation of the literature, for example Refs. [@Satow:2014lva; @Hong:2010at]. The Boltzmann equation reads $$\begin{aligned}
\begin{split} 
&v\cdot\partial_X n_{\pm f}({{\mathbf{k}}}, X) 
\pm eq_f\left({{\mathbf{E}}}+{{\mathbf{v}}}\times{{\mathbf{B}}}\right)(X) \cdot\nabla_{{{\mathbf{k}}}} n_{\pm f}({{\mathbf{k}}}, X)\\
&= C[n_{\pm f}], 
\end{split}\end{aligned}$$ where $n_{\pm f}({{\mathbf{k}}}, X)$ is the distribution function for the quark (antiquark) with momentum ${{\mathbf{k}}}$ at point $X$, and $v^\mu\equiv (1,{{\mathbf{v}}})$ with ${{\mathbf{v}}}\equiv {{\mathbf{k}}}/|{{\mathbf{k}}}|$. $C[n_{\pm f}]$ represents the collision effect among the quarks, which is given later.

Now we consider the situation in which the system at equilibrium is disturbed by weak external EM fields, so that the distribution function slightly deviates from the equilibrium, $n_{\pm f}({{\mathbf{k}}}, X)={n_F}(|{{\mathbf{k}}}|)+\delta n_{\pm f}({{\mathbf{k}}}, X)$. By linearizing the Boltzmann equation in terms of $\delta n_{\pm f}({{\mathbf{k}}}, X)$ and EM fields, we get $$\begin{aligned}
\label{eq:EOM-1}
\begin{split} 
v\cdot\partial_X \delta n_{\pm f}({{\mathbf{k}}}, X) 
\pm eq_f {{\mathbf{E}}}(X) \cdot{{\mathbf{v}}}{n_F}'(|{{\mathbf{k}}}|)
&= \delta C[n_{\pm f}], 
\end{split}\end{aligned}$$ where the magnetic field term disappears due to the isotropy of the system at equilibrium. $\delta C[n_{\pm f}]$ is a linearized form of $C[n_{\pm f}]$, whose expression reads $$\begin{aligned}
\delta C[n_{\pm f}]
&= -\tau^{-1} \left(\delta n_{\pm f}({{\mathbf{k}}}, X) 
\mp{n_F}'(|{{\mathbf{k}}}|) \delta\mu_f(X)
\right), \end{aligned}$$ in the relaxation time approximation. Here we have introduced the relaxation time $\tau\sim(g^4T)^{-1}$, and the shift of the chemical potential ($\delta\mu_f$) caused by the EM fields. The second term in the expression above is necessary, since the deviation of the distribution created by the shift of the chemical potential does not relax. The shift of the chemical potential is determined by the conservation law of particle number, $$\begin{aligned}
0&=\int\frac{d^3{{\mathbf{k}}}}{(2\pi)^3}\sum_{s=\pm 1} s  \delta C[n_{s f}] ,\end{aligned}$$ which reduces to $$\begin{aligned}
\label{eq:EOM-2}
\delta\mu_f(X)
&= -\frac{1}{\chi}
\int\frac{d^3{{\mathbf{k}}}}{(2\pi)^3}\sum_{s=\pm 1} s  \delta n_{s f}({{\mathbf{k}}}, X).\end{aligned}$$ Here, $\chi\equiv T^2/6$. From Eqs. (\[eq:EOM-1\]) and (\[eq:EOM-2\]), we obtain the solutions, $$\begin{aligned}
\delta\mu_f(p)
&= ieq_f\frac{{{\mathbf{E}}}(p)\cdot \hat{{{\mathbf{p}}}}}{|{{\mathbf{p}}}|} \frac{1-(\omega+i\tau^{-1})A(p)}{1-i\tau^{-1}A(p)}, \\
\label{eq:Boltzmanneq-solution}
\delta n_{s f}({{\mathbf{k}}}, p)
&= -is{n_F}'(|{{\mathbf{k}}}|)\frac{eq_f{{\mathbf{E}}}(p)\cdot{{\mathbf{v}}}-\tau^{-1}\delta\mu_f(p)}{v\cdot p+i\tau^{-1}},\end{aligned}$$ where we have performed the Fourier transformation $(X\rightarrow p)$ and introduced $A(p)\equiv \ln[(\omega+|{{\mathbf{p}}}|+i\tau^{-1})/(\omega-|{{\mathbf{p}}}|+i\tau^{-1})]/(2|{{\mathbf{p}}}|)$.

The induced current is given by $$\begin{aligned}
j^\mu(p)
&= 2e\sum_f q_f{N_c}\int\frac{d^3{{\mathbf{k}}}}{(2\pi)^3} v^\mu
\sum_{s=\pm 1} s\delta n_{sf}({{\mathbf{k}}},p).\end{aligned}$$ By using this expression in momentum space and the linear response theory relation of Eq.(\[eq:linear-response\]), Eq. (\[eq:Boltzmanneq-solution\]) leads to the retarded Green function $G^R_{\mu\nu}(p)$, as is shown below. We note that $j^0=-2e\sum_f q_f{N_c}\chi\delta\mu_f$, which indicates that $\chi$ is essentially the susceptibility.

Transverse channel
------------------

The transverse component of $G^R$ is given by $$\begin{aligned}
\begin{split}
G_T(p)
&=4 {C_{\text{em}}}{N_c}\frac{1}{(2\pi)^3}
\int^{2\pi}_0 d\phi \int^1_{-1} d\cos\theta \int^\infty_0 d|{{\mathbf{k}}}| |{{\mathbf{k}}}|^2 \\
&~~~\times{n_F}'(|{{\mathbf{k}}}|)\frac{\omega\sin^2\theta\cos^2\phi}{\omega-|{{\mathbf{p}}}|\cos\theta+i\tau^{-1}} \\
&= -{C_{\text{em}}}{N_c}\chi\frac{\omega}{{{\mathbf{p}}}^2} \\
&~~~\times \left[\omega+i\tau^{-1}-A(p)(p^2+2i\tau^{-1}\omega+(i\tau^{-1})^2)\right].
\end{split}
\label{eq:GTfull}\end{aligned}$$ In the hydro limit, $\omega,|{{\mathbf{p}}}|\ll\tau^{-1}$, this reduces to $$\begin{aligned}
\begin{split}
G_T(p)
&\simeq {C_{\text{em}}}{N_c}\chi \frac{2}{3}\tau \omega
 \left[i-\tau\omega\right],
\end{split} \end{aligned}$$ where we have used $$\begin{aligned}
\begin{split}
A(p) & \simeq \frac{1}{i\tau^{-1}}
\left[1-\frac{\omega}{i\tau^{-1}}+\frac{3\omega^2+{{\mathbf{p}}}^2}{3(i\tau^{-1})^2}
-\frac{\omega(\omega^2+{{\mathbf{p}}}^2)}{(i\tau^{-1})^3}
\right].
\end{split}\end{aligned}$$ By comparing this expression with the hydro result of Eq.(\[eq:hydro-GR-T\]), we derive the following expression for the transport coefficients: $$\begin{aligned}
\sigma &= {C_{\text{em}}}{N_c}\chi \frac{2}{3}\tau
= {C_{\text{em}}}{N_c}\frac{T^2}{9}\tau,\\
\tau_J&=\tau, \\
\kappa_B&=0.\end{aligned}$$

Next, we obtain the spectral function for small $|{{\mathbf{p}}}|$. Expanding Eq.(\[eq:GTfull\]) in terms of $|{{\mathbf{p}}}|$, we derive $$\begin{aligned}
\begin{split}
G_T(p)
&\simeq -{C_{\text{em}}}{N_c}\chi \frac{2}{3}\frac{\omega}{\omega+i\tau^{-1}} 
 \left[1+\frac{1}{5}\frac{{{\mathbf{p}}}^2}{(\omega+i\tau^{-1})^2}\right],
\end{split} 
\label{eq:GTexpanded}\end{aligned}$$ where we have used $$\begin{aligned}
\begin{split}
& A(p) \simeq \\
& \frac{1}{\omega+i\tau^{-1}}
\left[1+\frac{1}{3}\left(\frac{|{{\mathbf{p}}}|}{\omega+i\tau^{-1}}\right)^2
 +\frac{1}{5}\left(\frac{|{{\mathbf{p}}}|}{\omega+i\tau^{-1}}\right)^4\right].
\end{split}\end{aligned}$$ The imaginary part of Eq.(\[eq:GTexpanded\]) reads $$\begin{aligned}
\label{eq:transport-T}
\begin{split}
\rho_T(p)
&\simeq {C_{\text{em}}}{N_c}\chi \frac{2}{3}  \frac{\tau^{-1} \omega}{\omega^2+\tau^{-2}}
 \left[1 +\frac{{{\mathbf{p}}}^2}{5} \frac{(3\omega^2-\tau^{-2})}{(\omega^2+\tau^{-2})^2}\right].
\end{split} \end{aligned}$$

Longitudinal channel {#longitudinal-channel}
--------------------

The longitudinal component of $G^R$ is given by $$\begin{aligned}
\begin{split}
G^R_{00}(p)
&= 2{C_{\text{em}}}{N_c}\chi\frac{1-(\omega+i\tau^{-1})A(p)}{1-i\tau^{-1}A(p)}.
\end{split}
\label{eq:G00full}\end{aligned}$$

In the hydro limit, $\omega,|{{\mathbf{p}}}|\ll\tau^{-1}$, this reduces to $$\begin{aligned}
\label{eq:G00-hydrolim}
\begin{split}
G^R_{00}(p)
&\simeq i{{\mathbf{p}}}^2 \frac{2{C_{\text{em}}}{N_c}\chi \tau}{3}
\frac{1}
{(\omega+i\tau{{\mathbf{p}}}^2/3)},
\end{split} \end{aligned}$$ where we have retained only the leading order terms. Comparing this expression with the result of hydrodynamics (\[eq:hydro-GR-L\]), we get $$\begin{aligned}
D&= \frac{\tau}{3}, \end{aligned}$$ and also confirm the Einstein relation, $$\begin{aligned}
\sigma
&= 2{C_{\text{em}}}{N_c}\chi D,\end{aligned}$$ where the factor 2${N_c}$ originates from the spin and color degrees of freedom of the quark.

Now, we can obtain the spectral function for small $|{{\mathbf{p}}}|$. Expanding Eq.(\[eq:G00full\]) in terms of $|{{\mathbf{p}}}|$, we derive $$\begin{aligned}
\begin{split}
G^R_{00}(p)
&\simeq -\frac{2}{3}{C_{\text{em}}}{N_c}\chi \frac{|{{\mathbf{p}}}|^2}{\omega+i\tau^{-1}}
\frac{1}{\omega} \\
&~~~\times \left[1+\left(\frac{3}{5}+\frac{i\tau^{-1}}{3\omega}\right)
\left(\frac{|{{\mathbf{p}}}|}{\omega+i\tau^{-1}}\right)^2
\right].
\end{split} \end{aligned}$$ The imaginary part of the above expression reads $$\begin{aligned}
\label{eq:transport-00}
\begin{split}
\rho^{00}(p)
&\simeq 
\frac{2}{3}{C_{\text{em}}}{N_c}\chi \frac{{{\mathbf{p}}}^2}{\omega} \frac{\tau^{-1}}{\omega^2+\tau^{-2}} \\
&~~~\times \left[1
+{{\mathbf{p}}}^2 \frac{2}{5}\left(\tau^{-2}+\frac{11}{3}\omega^2\right)
\frac{1}{(\omega^2+\tau^{-2})^2}
\right].
\end{split} \end{aligned}$$

Evaluation of the continuum at weak coupling {#app:continuum}
============================================

In this appendix, we evaluate the continuum in the weak coupling and massless limit. In the free limit, the Green function of the EM current can be calculated by using Wick’s theorem as [@Altherr:1989jc] $$\begin{aligned}
\begin{split} 
G^R_{\mu\nu}(x)
&= i\theta(t) C_{\text{em}}{N_c}{\mathrm{Tr}}\Bigl[\gamma_\mu S^>(x)\gamma_\nu S^<(-x) \\
&~~~-\gamma_\mu S^<(x)\gamma_\nu S^>(-x)
\Bigr],
\end{split}\end{aligned}$$ where $S^>(x)\equiv \langle \psi(x)\overline{\psi}(0) \rangle$ and $S^<(x)\equiv \langle \overline{\psi}(0) \psi(x)\rangle$. Here, we have omitted the flavor indices for simplicity. By performing the Fourier transformation and taking the imaginary part, the spectral function reads $$\begin{aligned}
\begin{split} 
\rho_{\mu\nu}(p)
&= \frac{1}{2} C_{\text{em}}{N_c}\int \frac{d^4k}{(2\pi)^4}
{\mathrm{Tr}}\Bigl[\gamma_\mu S^>(k)\gamma_\nu S^<(k-p) \\
&~~~-\gamma_\mu S^<(k)\gamma_\nu S^>(k-p)
\Bigr]\\
&=2 C_{\text{em}}{N_c}\int \frac{d^4k}{(2\pi)^4}
{\mathrm{Tr}}\Bigl[\gamma_\mu {\ooalign{\hfil/\hfil\crcr$k$}}\gamma_\nu ({\ooalign{\hfil/\hfil\crcr$k$}}-{\ooalign{\hfil/\hfil\crcr$p$}}) \Bigr] \\
&~~~\times \rho^0(k) \rho^0(k-p)
[{n_F}(k^0-\omega)-{n_F}(k^0)],
\end{split}
\label{eq:rhomunu1}\end{aligned}$$ where we have used $S^>(k)= {\ooalign{\hfil/\hfil\crcr$k$}}2\rho^0(k)[1-{n_F}(k^0)]$ and $S^<(k)={\ooalign{\hfil/\hfil\crcr$k$}}2\rho^0(k){n_F}(k^0)$, and introduced the free quark spectral function as $\rho^0(k)\equiv \pi{\mathrm{sgn}}(k^0)\delta(k^2)$. The two delta functions can be written as $$\begin{aligned}
\label{eq:deltafunc-1}
\delta(k^2)&=\sum_{s=\pm 1}\frac{\delta(k^0-s|{{\mathbf{k}}}|)}{2|{{\mathbf{k}}}|},\\
\nonumber
\delta([k-p]^2)&= \frac{1}{2|{{\mathbf{k}}}||{{\mathbf{p}}}|}\delta\left[\cos\theta-\frac{2k^0\omega-p^2}{2|{{\mathbf{k}}}||{{\mathbf{p}}}|}\right] \\
\label{eq:deltafunc-2}
&~~~\times\theta\left[-p^2\left\{(k^0)^2-k^0\omega+\frac{p^2}{4}\right\}\right],\end{aligned}$$ where we adopted the standard polar coordinate, by assigning ${{\mathbf{p}}}$ to point into the $z$-direction. We also note that $\theta$ in the first line is the angle between ${{\mathbf{k}}}$ and ${{\mathbf{p}}}$ while that in the second line is a step function. Therefore, the expression of Eq.(\[eq:rhomunu1\]) becomes $$\begin{aligned}
\begin{split}
\rho_{\mu\nu}(p)
&= C_{\text{em}}{N_c}\frac{1}{8\pi^2|{{\mathbf{p}}}|}
\int^{2\pi}_0 d\phi \int^\infty_0 d|{{\mathbf{k}}}| \sum_{s=\pm 1} \\
&~~~\times \Bigl[k_\mu (k-p)_\nu+k_\nu (k-p)_\mu+g_{\mu\nu} k\cdot p \Bigr] \\
&~~~\times  {\mathrm{sgn}}(s|{{\mathbf{k}}}|-\omega)
s \theta\left[-\left\{(k^0)^2-k^0\omega+\frac{p^2}{4}\right\}\right]\\
&~~~\times [{n_F}(s|{{\mathbf{k}}}|-\omega)-{n_F}(s|{{\mathbf{k}}}|)],
\end{split}
\label{eq:rhomunu2}\end{aligned}$$ where the values for $k^0$ and $\cos\theta$ are determined by the delta functions of Eqs. (\[eq:deltafunc-1\]) and (\[eq:deltafunc-2\]), and we have used $p^2>0$, which is justified because we consider the case that $\omega \sim T\gg |{{\mathbf{p}}}|$. From now on, we focus on the case for which $\omega>0$. Then, only the contribution with $s=+1$ remains, and the step function restricts the range of $|{{\mathbf{k}}}|$ to $k_-<|{{\mathbf{k}}}|<k_+$, where $k_\pm\equiv (\omega\pm|{{\mathbf{p}}}|)/2$. Thus, Eq.(\[eq:rhomunu2\]) becomes $$\begin{aligned}
\label{eq:rho-continuum}
\begin{split}
\rho_{\mu\nu}(p)
&= -C_{\text{em}}{N_c}\frac{1}{8\pi^2|{{\mathbf{p}}}|}
\int^{2\pi}_0 d\phi \int^{k_+}_{k_-} d|{{\mathbf{k}}}|  \\
&~~~\times \Bigl[k_\mu (k-p)_\nu+k_\nu (k-p)_\mu+g_{\mu\nu} k\cdot p \Bigr] \\
&~~~\times   [1-{n_F}(-|{{\mathbf{k}}}|+\omega)-{n_F}(|{{\mathbf{k}}}|)],
\end{split}\end{aligned}$$ where we have used $\omega-|{{\mathbf{k}}}|>(\omega-|{{\mathbf{p}}}|)/2>0$. From the distribution factor $[1-{n_F}(-|{{\mathbf{k}}}|+\omega)-{n_F}(|{{\mathbf{k}}}|)]=[1-{n_F}(-|{{\mathbf{k}}}|+\omega)][1-{n_F}(|{{\mathbf{k}}}|)]-{n_F}(-|{{\mathbf{k}}}|+\omega){n_F}(|{{\mathbf{k}}}|)$, we see that the physical processes corresponding to this expression are the quark antiquark pair-creation process and its inverse.

Transverse channel
------------------

$\delta\rho^T(p)$ can be evaluated by setting $\mu=\nu=1$ and subtracting the $T=0$ part from Eq. (\[eq:rho-continuum\]): $$\begin{aligned}
\begin{split}
\delta\rho_{T}(p)
&= C_{\text{em}}{N_c}\frac{1}{8\pi^2|{{\mathbf{p}}}|}
\int^{2\pi}_0 d\phi \int^{k_+}_{k_-} d|{{\mathbf{k}}}|  \\
&~~~\times \Bigl[2{{\mathbf{k}}}^2\cos^2\phi(1-\cos^2\theta)-|{{\mathbf{k}}}|\omega+|{{\mathbf{k}}}||{{\mathbf{p}}}|\cos\theta \Bigr] \\
&~~~\times   [{n_F}(-|{{\mathbf{k}}}|+\omega)+{n_F}(|{{\mathbf{k}}}|)]\\
&= -C_{\text{em}}{N_c}\frac{p^2}{4\pi|{{\mathbf{p}}}|^3}
\int^{|{{\mathbf{p}}}|/2}_{-|{{\mathbf{p}}}|/2} d|{{\mathbf{k}}}|  
\left(\frac{{{\mathbf{p}}}^2}{4}+{{\mathbf{k}}}^2\right) \\
&~~~\times \left[{n_F}\left(\frac{\omega}{2}-|{{\mathbf{k}}}|\right)+{n_F}\left(\frac{\omega}{2}+|{{\mathbf{k}}}|\right)\right],
\end{split}\end{aligned}$$ where in the last line, we have used $\cos\theta=(2|{{\mathbf{k}}}|\omega-p^2)(2|{{\mathbf{k}}}||{{\mathbf{p}}}|)$ and have changed the integration variable as $|{{\mathbf{k}}}|\rightarrow |{{\mathbf{k}}}|-\omega/2$.

We can safely expand this in terms of $|{{\mathbf{k}}}|/\omega$, because $|{{\mathbf{k}}}|\simeq |{{\mathbf{p}}}|\ll \omega$, which leads to $$\begin{aligned}
\label{eq:continuum-T}
\begin{split}
& \delta\rho_{T}(p) \\ 
&\simeq - C_{\text{em}}{N_c}\frac{\omega^2}{\pi|{{\mathbf{p}}}|^3}\left(1-\frac{{{\mathbf{p}}}^2}{\omega^2}\right)
\int^{|{{\mathbf{p}}}|/2}_{0} d|{{\mathbf{k}}}|  
\left(\frac{{{\mathbf{p}}}^2}{4}+{{\mathbf{k}}}^2\right) \\
&~~~\times \left[{n_F}\left(\frac{\omega}{2}\right)+\frac{{{\mathbf{k}}}^2}{2}{n_F}''\left(\frac{\omega}{2}\right)\right]\\
&\simeq - C_{\text{em}}{N_c}\frac{\omega^2}{6\pi}\left(1-\frac{{{\mathbf{p}}}^2}{\omega^2}\right) 
 \left[{n_F}\left(\frac{\omega}{2}\right)+\frac{{{\mathbf{p}}}^2}{20}{n_F}''\left(\frac{\omega}{2}\right)\right]\\
 &\simeq -C_{\text{em}}{N_c}\frac{\omega^2}{6\pi}
 \Biggl[{n_F}\left(\frac{\omega}{2}\right) \\
&~~~ +{{\mathbf{p}}}^2\left\{\frac{1}{20}{n_F}''\left(\frac{\omega}{2}\right)
 -\frac{1}{\omega^2}{n_F}\left(\frac{\omega}{2}\right)\right\}\Biggr] .
\end{split}\end{aligned}$$

Longitudinal channel {#longitudinal-channel-1}
--------------------

By setting $\mu=\nu=0$ in Eq. (\[eq:rho-continuum\]), we get $$\begin{aligned}
\begin{split}
\delta\rho_{00}(p)
&= C_{\text{em}}{N_c}\frac{1}{2\pi|{{\mathbf{p}}}|}
\int^{|{{\mathbf{p}}}|/2}_{-|{{\mathbf{p}}}|/2} d|{{\mathbf{k}}}|  
\left[ {{\mathbf{k}}}^2- \frac{{{\mathbf{p}}}^2}{4} \right] \\
&~~~\times   \left[{n_F}\left(\frac{\omega}{2}-|{{\mathbf{k}}}|\right)+{n_F}\left(\frac{\omega}{2}+|{{\mathbf{k}}}|\right)\right],
\end{split}\end{aligned}$$ where we have changed the integration variable as before.

Expanding the integrand in terms of $|{{\mathbf{p}}}|/\omega$, we derive $$\begin{aligned}
\label{eq:continuum-00}
\begin{split}
\delta\rho_{00}(p)
&= C_{\text{em}}{N_c}\frac{2}{\pi|{{\mathbf{p}}}|}
\int^{|{{\mathbf{p}}}|/2}_{0} d|{{\mathbf{k}}}|  
\left[ {{\mathbf{k}}}^2- \frac{{{\mathbf{p}}}^2}{4} \right] \\
&~~~\times \left[{n_F}\left(\frac{\omega}{2}\right)+\frac{{{\mathbf{k}}}^2}{2}{n_F}''\left(\frac{\omega}{2}\right)\right]\\
&= -C_{\text{em}}{N_c}\frac{{{\mathbf{p}}}^2}{6\pi} 
 \left[{n_F}\left(\frac{\omega}{2}\right)+\frac{{{\mathbf{p}}}^2}{40}{n_F}''\left(\frac{\omega}{2}\right)\right] .
\end{split}  \end{aligned}$$

Evaluation of the UV tail at weak coupling {#app:UVtail}
==========================================

In this appendix, we briefly recapitulate the derivation of the UV tail in the EM current spectral function from the OPE [@CaronHuot:2009ns]. The UV behavior of the EM current retarded correlator is described by the OPE of Eqs.(\[eq:OPE-wo-RG-T\]) and (\[eq:OPE-wo-RG-L\]). Among the three terms in these expressions, only $\langle T^{00}_f\rangle$ is not RG invariant. This operator yields imaginary parts of the retarded correlator, as can be understood as follows: The scaling relation (\[eq:scaling-T’-Ttilde\]) can be rewritten as $$\begin{aligned}
\begin{split} 
T'{}^{00}_{f}(\kappa)
&\simeq  T'{}^{00}_{f}(\kappa_0)
+a'\ln\left(\frac{\kappa^2_0}{\kappa^2}\right)
\frac{b_0}{4\pi}\alpha_sT'{}^{00}_{f} ,\\
 \tilde{T}^{00}(\kappa)
&\simeq  \tilde{T}^{00}(\kappa_0)
+\tilde{a}\ln\left(\frac{\kappa^2_0}{\kappa^2}\right)
\frac{b_0}{4\pi}\alpha_s\tilde{T}^{00}_{f},
\end{split} \end{aligned}$$ when $\kappa$ is close to $\kappa_0$. It was shown in Ref.[@CaronHuot:2009ns] that the factor $\ln(\kappa^2_0/\kappa^2)$ generates an imaginary contribution $i\pi$, due to the analytic continuation to the real time. Following this prescription, the imaginary parts of the retarded correlators (\[eq:OPE-wo-RG-T\]) and (\[eq:OPE-wo-RG-L\]) read $$\begin{aligned}
\nonumber
\delta \rho_T(p) 
&=e^2\sum q^2_f \frac{8}{9}\frac{\omega^2+{{\mathbf{p}}}^2}{(p^2)^2}
\alpha_s(\omega) \\
\label{eq:deltarho_T_UVtail}
&~~~\times  \left( 2{C_F}\delta\left\langle T'{}^{00}_{f}(\omega)  \right\rangle
+\frac{1}{{N_f}}  \delta\left\langle \tilde{T}^{00}(\omega) \right\rangle \right) ,\\
\nonumber
\delta \rho_{00}(p)
&=  e^2\sum q^2_f \frac{8}{9}\frac{{{\mathbf{p}}}^2}{(p^2)^2}
\alpha_s(\omega) \\
&~~~\times  \left( 2{C_F}\delta\left\langle T'{}^{00}_{f}(\omega)  \right\rangle
+\frac{1}{{N_f}}  \delta\left\langle \tilde{T}^{00}(\omega) \right\rangle \right) .\end{aligned}$$ We note that this expression is valid when the OPE is reliable ($\omega\gg T,\Lambda_{\text {QCD}}$).

In the chiral and weak coupling limits, the operator expectation values at the renormalization scale $\kappa_0\sim T$ read $$\begin{aligned}
\label{eq:T00f-free}
\langle T^{00}_f\rangle&= {N_c}\frac{7\pi^2T^4}{60},\\
\label{eq:T00g-free}
\langle T^{00}_g\rangle&=2{C_F}{N_c}\frac{\pi^2T^4}{15},\end{aligned}$$ which, by using the scaling relation of Eq.(\[eq:scaling-T’-Ttilde\]), leads to $$\begin{aligned}
\nonumber
\delta \rho_T(p) 
&={C_{\text{em}}}\frac{1}{\omega^2} 
\left(1+3\frac{{{\mathbf{p}}}^2}{\omega^2}\right)
\alpha_s(\kappa_0) {N_c}{C_F}\frac{4\pi^2T^4}{27} \\
\label{eq:UVtail-T}
&~~~\times  \left[\frac{\ln\left(\kappa_0/\Lambda_{\text{QCD}} \right)}{\ln\left(\omega/\Lambda_{\text{QCD}} \right)}\right]^{\tilde{a}+1}  ,\\
\nonumber
\delta \rho_{00}(p)
&={C_{\text{em}}}\frac{{{\mathbf{p}}}^2}{\omega^4} 
\left(1+2\frac{{{\mathbf{p}}}^2}{\omega^2}\right)
\alpha_s(\kappa_0) {N_c}{C_F}\frac{4\pi^2T^4}{27}\\
\label{eq:UVtail-00}
&~~~\times  \left[\frac{\ln\left(\kappa_0/\Lambda_{\text{QCD}} \right)}{\ln\left(\omega/\Lambda_{\text{QCD}} \right)}\right]^{\tilde{a}+1}. \end{aligned}$$ Here, we have retained terms up to next-to-leading order in the small $|{{\mathbf{p}}}|$ expansion.

[99]{}

A. Adare [*et al.*]{} \[PHENIX Collaboration\], Phys. Rev. C [**91**]{}, no. 6, 064904 (2015) \[arXiv:1405.3940 \[nucl-ex\]\]; A. Adare [*et al.*]{} \[PHENIX Collaboration\], Phys. Rev. C [**94**]{}, no. 6, 064901 (2016) \[arXiv:1509.07758 \[nucl-ex\]\]; J. Adam [*et al.*]{} \[ALICE Collaboration\], Phys. Lett. B [**754**]{}, 235 (2016) \[arXiv:1509.07324 \[nucl-ex\]\]; L. Adamczyk [*et al.*]{} \[STAR Collaboration\], arXiv:1607.01447 \[nucl-ex\].

X. G. Huang, Rept. Prog. Phys.  [**79**]{}, no. 7, 076302 (2016) \[arXiv:1509.04073 \[nucl-th\]\]; K. Hattori and X. G. Huang, arXiv:1609.00747 \[nucl-th\]. L. McLerran and V. Skokov, Nucl. Phys. A [**929**]{}, 184 (2014) \[arXiv:1305.0774 \[hep-ph\]\]. K. Tuchin, Phys. Rev. C [**88**]{}, no. 2, 024911 (2013) \[arXiv:1305.5806 \[hep-ph\]\]. K. Tuchin, Phys. Rev. C [**93**]{}, no. 1, 014905 (2016) \[arXiv:1508.06925 \[hep-ph\]\].

V. Bernard and U. G. Meissner, Nucl. Phys. A [**489**]{}, 647 (1988); G. E. Brown and M. Rho, Phys. Rev. Lett.  [**66**]{}, 2720 (1991); T. Hatsuda and S. H. Lee, Phys. Rev. C [**46**]{}, no. 1, R34 (1992).

R. Baier, B. Pire and D. Schiff, Phys. Rev. D [**38**]{}, 2814 (1988); M. Laine, JHEP [**1311**]{}, 120 (2013) \[arXiv:1310.0164 \[hep-ph\]\]. S. Caron-Huot, P. Kovtun, G. D. Moore, A. Starinets and L. G. Yaffe, JHEP [**0612**]{}, 015 (2006) \[hep-th/0607237\]. C. Gale, Y. Hidaka, S. Jeon, S. Lin, J.-F. Paquet, R. D. Pisarski, D. Satow and V. V. Skokov [*et al.*]{}, Phys. Rev. Lett.  [**114**]{}, 072301 (2015) \[arXiv:1409.4778 \[hep-ph\]\]; Y. Hidaka, S. Lin, R. D. Pisarski and D. Satow, JHEP [**1510**]{}, 005 (2015) \[arXiv:1504.01770 \[hep-ph\]\]; D. Satow and W. Weise, Phys. Rev. D [**92**]{}, 056001 (2015) \[arXiv:1505.03869 \[hep-ph\]\]; C. A. Islam, S. Majumder, N. Haque and M. G. Mustafa, JHEP [**1502**]{}, 011 (2015) \[arXiv:1411.6407 \[hep-ph\]\]. P. Gubler and W. Weise, Phys. Lett. B [**751**]{}, 396 (2015) \[arXiv:1507.03769 \[hep-ph\]\]; P. Gubler and W. Weise, Nucl. Phys. A [**954**]{}, 125 (2016) \[arXiv:1602.09126 \[hep-ph\]\]; T. Hayata, Phys. Rev. D [**88**]{}, no. 3, 036007 (2013) \[arXiv:1306.3759 \[hep-ph\]\]; P. M. Hohler and R. Rapp, Nucl. Phys. A [**892**]{}, 58 (2012) \[arXiv:1204.6309 \[hep-ph\]\]; Phys. Lett. B [**731**]{}, 103 (2014) \[arXiv:1311.2921 \[hep-ph\]\]; P. Gubler and K. Ohtani, Phys. Rev. D [**90**]{}, 094002 (2014) \[arXiv:1404.7701 \[hep-ph\]\]; P. Gubler, K. Morita and M. Oka, Phys. Rev. Lett. [**107**]{}, 092003 (2011) \[arXiv:1104.4436 \[hep-ph\]\]; K. Suzuki, P. Gubler, K. Morita and M. Oka, Nucl. Phys. A [**897**]{}, 28 (2013) \[arXiv:1204.1173 \[hep-ph\]\]; C. A. Dominguez, M. Loewe, J. C. Rojas and Y. Zhang, Phys. Rev. D [**81**]{}, 014007 (2010) \[arXiv:0908.2709 \[hep-ph\]\]; K. Morita and S. H. Lee, Phys. Rev. Lett. [**100**]{}, 022301 (2008) \[arXiv:0704.2021 \[nucl-th\]\]; A. Ayala, C. A. Dominguez and M. Loewe, Adv. High Energy Phys.  [**2017**]{}, 9291623 (2017) \[arXiv:1608.04284 \[hep-ph\]\]. J. I. Kapusta and E. V. Shuryak, Phys. Rev. D [**49**]{}, 4694 (1994) \[hep-ph/9312245\]. S. Z. Huang and M. Lissia, Phys. Lett. B [**348**]{}, 571 (1995) \[hep-ph/9404275\]; Phys. Rev. D [**52**]{}, 1134 (1995) \[hep-ph/9412246\].

S. Zschocke, O. P. Pavlenko and B. Kampfer, Eur. Phys. J. A [**15**]{}, 529 (2002) \[nucl-th/0205057\].

H. Kim, P. Gubler and S. H. Lee, Phys. Lett. B [**772**]{}, 194 (2017) \[arXiv:1703.04848 \[hep-ph\]\]. G. Chanfray, R. Rapp and J. Wambach, Phys. Rev. Lett. [**76**]{}, 368 (1996) \[hep-ph/9508353\].

F. Klingl, N. Kaiser and W. Weise, Nucl. Phys. A [**624**]{}, 527 (1997) \[hep-ph/9704398\].

H. -T. Ding, A. Francis, O. Kaczmarek, F. Karsch, E. Laermann and W. Soeldner, Phys. Rev. D [**83**]{}, 034504 (2011) \[arXiv:1012.4963 \[hep-lat\]\]. A. Amato, G. Aarts, C. Allton, P. Giudice, S. Hands and J. -I. Skullerud, Phys. Rev. Lett.  [**111**]{}, 172001 (2013) \[arXiv:1307.6763 \[hep-lat\]\]. G. Aarts, C. Allton, J. Foley, S. Hands and S. Kim, Phys. Rev. Lett.  [**99**]{}, 022002 (2007) \[hep-lat/0703008 \[HEP-LAT\]\]. B. B. Brandt, A. Francis, B. J$\mathrm{\ddot{a}}$ger and H. B. Meyer, Phys. Rev. D [**93**]{}, no. 5, 054510 (2016) \[arXiv:1512.07249 \[hep-lat\]\]. H. B. Meyer, Eur. Phys. J. A [**47**]{}, 86 (2011) \[arXiv:1104.3708 \[hep-lat\]\]. M. Jarrell and J. E. Gubernatis, Phys. Rept.  [**269**]{}, 133 (1996). M. Asakawa, T. Hatsuda and Y. Nakahara, Prog. Part. Nucl. Phys.  [**46**]{}, 459 (2001) \[hep-lat/0011040\]. P. Gubler and M. Oka, Prog. Theor. Phys.  [**124**]{}, 995 (2010) \[arXiv:1005.2459 \[hep-ph\]\].

G. Backus and F. Gilbert, Geophys. J. R. Astron. Soc. [**16**]{}, 169 (1968).

G. Backus and F. Gilbert, Phil. Trans. R. Soc. A [**266**]{}, 123 (1970).

B. B. Brandt, A. Francis, H. B. Meyer and D. Robaina, Phys. Rev. D [**92**]{}, no. 9, 094510 (2015) \[arXiv:1506.05732 \[hep-lat\]\]. L. Schlessinger, Phys. Rev. [**167**]{}, 1411 (1968).

R. A. Tripolt, I. Haritan, J. Wambach and N. Moiseyev, arXiv:1610.03252 \[hep-ph\]. R. A. Tripolt, P. Gubler, M. Ulybyshev and L. von Smekal, in preparation.

P. Gubler and D. Satow, Phys. Rev. D [**94**]{}, no. 9, 094042 (2016) \[arXiv:1602.08265 \[hep-ph\]\]. P. Romatschke and D. T. Son, Phys. Rev. D [**80**]{}, 065021 (2009) \[arXiv:0903.3946 \[hep-ph\]\]. A. Ayala, C. A. Dominguez, L. A. Hern$\mathrm{\acute{a}}$ndez, M. Loewe and A. J. Mizher, Phys. Rev. D [**88**]{}, 114028 (2013) \[arXiv:1309.4135 \[hep-ph\]\]. S. Caron-Huot, Phys. Rev. D [**79**]{}, 125009 (2009) \[arXiv:0903.3958 \[hep-ph\]\].

M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B [**147**]{}, 448 (1979). M. E. Peskin and D. V. Schroeder, Reading, USA: Addison-Wesley (1995) 842 p.

P. Kovtun, J. Phys. A [**45**]{}, 473001 (2012) \[arXiv:1205.5040 \[hep-th\]\]. P. Kovtun and L. G. Yaffe, Phys. Rev. D [**68**]{}, 025007 (2003) \[hep-th/0303010\].

B. B. Brandt, A. Francis and H. B. Meyer, Phys. Rev. D [**89**]{}, no. 3, 034506 (2014) \[arXiv:1310.5160 \[hep-lat\]\]. D. Bernecker and H. B. Meyer, Eur. Phys. J. A [**47**]{}, 148 (2011) \[arXiv:1107.4388 \[hep-lat\]\].

G. D. Moore and J. -M. Robert, hep-ph/0607172. E. Braaten, R. D. Pisarski and T. -C. Yuan, Phys. Rev. Lett.  [**64**]{}, 2242 (1990).

A. I. Bochkarev and M. E. Shaposhnikov, Nucl. Phys. B [**268**]{}, 220 (1986). A. Bazavov [*et al.*]{} \[HotQCD Collaboration\], Phys. Rev. D [**90**]{}, 094503 (2014) \[arXiv:1407.6387 \[hep-lat\]\].

A. Francis and O. Kaczmarek, Prog. Part. Nucl. Phys.  [**67**]{}, 212 (2012) \[arXiv:1112.4802 \[hep-lat\]\]. O. Kaczmarek, E. Laermann, M. M$\mathrm{\ddot{u}}$ller, F. Karsch, H. T. Ding, S. Mukherjee, A. Francis and W. Soeldner, PoS ConfinementX [****]{}, 185 (2012) \[arXiv:1301.7436 \[hep-lat\]\]. H. T. Ding, O. Kaczmarek and F. Meyer, PoS LATTICE [**2014**]{}, 216 (2015) \[arXiv:1412.5869 \[hep-lat\]\]. B. B. Brandt, A. Francis, H. B. Meyer and H. Wittig, JHEP [**1303**]{}, 100 (2013) \[arXiv:1212.4200 \[hep-lat\]\]. Y. Kwon, M. Procura and W. Weise, Phys. Rev. C [**78**]{}, 055203 (2008) \[arXiv:0803.3262 \[nucl-th\]\]. J. Hong and D. Teaney, Phys. Rev. C [**82**]{}, 044908 (2010) \[arXiv:1003.0699 \[nucl-th\]\]. D. Satow and H. U. Yee, Phys. Rev. D [**90**]{}, 014027 (2014) \[arXiv:1406.1150 \[hep-ph\]\]. T. Altherr and P. Aurenche, Z. Phys. C [**45**]{}, 99 (1989).

[^1]: This contour actually runs slightly above the real axis, so that it does not overlap with the singularities such as the continuum or the diffusion pole at zero momentum, which appear on it.

[^2]: We note that our convention for $\rho_{\mu\nu}$ differs from the popular $\rho_{\mu\nu}=2{\text{Im}}G^R_{\mu\nu}$ by a factor of 2.

[^3]: We could also choose $\kappa^2=p^2$, which however would not change the results of this paper.
