\section{Introduction and statement of results}
\label{sec:intro}}

The \emph{Riemann zeta function} is a central
object of study in analytic number theory. In terms of
the complex parameter $s=\sigma+it$, the zeta function is defined in the
half-plane $\sigma>1$ by two equivalent expressions:
$$
\zeta(s)\defeq\sum_{n\ge 1} n^{-s}
=\prod_{p\text{~prime}}(1-p^{-s})^{-1}.
$$
Riemann~\cite{Riemann} showed
that $\zeta(s)$ extends analytically to a meromorphic
function in the whole complex plane, its only singularity being a
simple pole at $s=1$. Moreover, the zeta function
satisfies a functional equation relating its values at
$s$ and $1-s$.\footnote{There are many excellent accounts of the
theory of the Riemann zeta function; we refer the reader to
Titchmarsh~\cite{Titchmarsh} and to Borwein \emph{et al}~\cite{Borwein}
for essential background.}
The \emph{Riemann Hypothesis} (RH) asserts that if $\rho=\beta+i\gamma$
is a zero of $\zeta(s)$ with real part $\beta>0$, then $\beta=\frac12$.

More generally, for a Dirichlet character $\chi\bmod q$,
the Dirichlet \text{$L$-function} $L(s,\chi)$ is defined
for $\sigma>1$ by:
$$
L(s,\chi)\defeq\sum_{n\ge 1} \chi(n)n^{-s}
=\prod_{p\text{~prime}}(1-\chi(p)p^{-s})^{-1}.
$$
The function $L(s,\chi)$ extends to a meromorphic
function (which is entire if $\chi$ is nonprincipal),
and when $\chi$ is primitive
it satisfies a simple functional equation relating its values at
$s$ and $1-s$; see, e.g., Bump~\cite[Chapter~1]{Bump}.
The \emph{Generalized Riemann Hypothesis} (GRH), which was perhaps first
formulated by Piltz in 1884 (see Davenport~\cite{Daven}), asserts that
if $\rho=\beta+i\gamma$ is a zero of $L(s,\chi)$ with
$\beta>0$, then $\beta=\frac12$.

For any \emph{principal} character $\chi_0\bmod q$ one has
$$
L(s,\chi_0)=\zeta(s)\prod_{p\,\mid\,q}(1-p^{-s}),
$$
hence RH is equivalent to GRH for $L(s,\chi_0)$. On the other hand,
for \emph{nonprincipal} characters $\chi$, no direct relationship between
RH and GRH for $L(s,\chi)$ has been previously established. To establish
such a connection, we study the following weak form of the GRH
for Dirichlet $L$-functions.

\bigskip\noindent{\sc Hypothesis ${\rm GRH}[\beta_0]$}:
{\it Given $\beta_0\in[\frac12,1)$, the
inequality $\beta\le\beta_0$ holds for all zeros $\rho=\beta+i\gamma$ of
an arbitrary Dirichlet $L$-function $L(s,\chi)$.}

\bigskip\noindent
Note that ${\rm GRH}[\tfrac12]$ is equivalent to that assertion that GRH
holds for all Dirichlet \text{$L$-functions}.
In the present paper, we show that hypothesis ${\rm GRH}[\frac{9}{10}]$
can be reformulated entirely in terms of certain distributional
properties of the zeros of the Riemann zeta function.

To state our results, we introduce some notation. If $T>0$ is not the
ordinate of a zero of the zeta function, then $N(T)$ is used to denote the number
of zeros $\rho=\beta+i\gamma$ of $\zeta(s)$ in the rectangle $0<\beta<1$,
$0<\gamma<T$, and we define
$$
S(T)\defeq\frac1\pi{\rm arg}\,\zeta(1/2+iT).
$$
If $\gamma>0$ is the ordinate of a zero, then we set
$$
N(\gamma)\defeq\frac12\{N(\gamma^+)+N(\gamma^-)\},\qquad
S(\gamma)\defeq\frac12\{S(\gamma^+)+S(\gamma^-)\}.
$$
Using an explicit form of the well known relation
(see, e.g., Montgomery and Vaughan\cite[Corollary 14.2]{MontVau})
\begin{equation}\label{eq:NTSTrelation}
N(T)=\frac{T}{2\pi}\log\frac{T}{2\pi}-\frac{T}{2\pi}+\frac{7}{8}
+S(T)+O(T^{-1})\qquad(T>0),
\end{equation}
one sees that the difference $S(\gamma^+)-S(\gamma^-)$ is an integer
for every zero $\rho=\beta+i\gamma$ of $\zeta(s)$ with $\gamma>0$. Extending
the definition of $S(T)$ appropriately, this also holds for complex zeros
with $\gamma<0$. Therefore, writing
$$
\e(u)\defeq \mathrm{e}^{2\pi iu}\qquad(u\in\R),
$$
we can unambiguously define
\begin{equation}\label{eq:Zrho-defn}
\EuScript{Z}(\rho)\defeq\lim\limits_{T\to\gamma}\overline{\e\big(S(T)\big)}
\end{equation}
(the common value of the left and right limits)
for every complex zero $\rho$ of $\zeta(s)$.

Our main results are as follows.
Throughout the paper, we use $C_c^\infty(\R^+)$ to denote the space of smooth
functions $f:\R^+\to\C$ with compact support in $\R^+$.

\bigskip

\begin{theorem}\label{thm:RHvsGRH}
Assume RH. Suppose that for every function $\EuScript{B}\in C_c^\infty(\R^+)$
and every rational number $\xi\in\Q\cap(0,1)$, the bound
\begin{equation}\label{eq:superbound}
\sum_{\rho=\frac12+i\gamma}\xi^{-i\gamma}\,\EuScript{Z}(\rho)
\,\EuScript{B}\Big(\frac{\gamma}{2\pi X}\Big)
\,\mathop{\ll}_{\xi,\EuScript{B},\eps}\,X^{9/10+\eps}
\end{equation}
holds for any $\eps>0$, where the sum runs over all
complex zeros $\rho=\frac12+i\gamma$ of $\zeta(s)$,
the quantity $\EuScript{Z}(\rho)$ is given by \eqref{eq:Zrho-defn},
and the implied constant depends only on
$\xi$, $\EuScript{B}$ and $\eps$. Then, the hypothesis ${\rm GRH}[\frac{9}{10}]$ is true.
\end{theorem}

We reiterate that the hypotheses of Theorem~\ref{thm:RHvsGRH}
involve only properties of the zeros of the Riemann zeta function.
\emph{No conditions are imposed on the zeros
of nonprincipal Dirichlet $L$-functions}.

We also have the following converse of Theorem~\ref{thm:RHvsGRH}.

\begin{theorem}\label{thm:RHvsGRH2}
Assume RH. If ${\rm GRH}[\frac{9}{10}]$ is true, then
for every $\EuScript{B}\in C_c^\infty(\R^+)$ and
every $\xi\in\Q\cap(0,1)$, the bound \eqref{eq:superbound}
holds for any $\eps>0$ with an implied constant that depends only on
$\xi$, $\EuScript{B}$ and $\eps$.
\end{theorem}

In view of Theorem~\ref{thm:RHvsGRH2}, it is natural to expect that the bound
\eqref{eq:superbound} holds for all $\EuScript{B}\in C_c^\infty(\R^+)$ and
$\xi\in\Q\cap(0,1)$; we conjecture that this is the case.

{\large\section{Approach}}

Using the explicit formula and assuming RH, we show that a sum of the form
\begin{equation}\label{eq:einstein}
\sum_{n\ge 1}\Lambda(n)\,\e(-n\xi)\,\EuScript{B}(n/X)
\end{equation}
with $\xi>0$ and $\EuScript{B}\in C_c^\infty(\R^+)$ is
equal to
\begin{equation}\label{eq:feynman}
-\ssum{\rho=\frac12+i\gamma}\xi^{-1/2-i\gamma}
\EuScript{Z}(\rho)\EuScript{B}\Big(\frac{\gamma}{2\pi\xi X}\Big)
\end{equation}
up to an error of at most $O_{\xi,\EuScript{B}}(X^{9/10})$;
see Theorem~\ref{thm:vonMangoldt-twist}.
To handle the integrals that
arise in our use of the explicit formula, we apply the method of the
stationary phase (however, some care is needed to
obtain adequate and explicit estimates for the error terms);
see Lemma~\ref{lem:the-child}, which is the main workhorse for
the proof of Theorem~\ref{thm:vonMangoldt-twist}.

Given the close relationship (under RH) between the sums \eqref{eq:einstein}
and \eqref{eq:feynman}, the bound \eqref{eq:superbound} for
rational numbers $\xi\in\Q\cap(0,1)$ allows us to deduce an equally strong bound
\begin{equation}\label{eq:thorne}
\sum_{n\ge 1}\Lambda(n)\chi(n)\EuScript{B}(n/X)
\,\mathop{\ll}_{q,\EuScript{B},\eps}\,X^{9/10+\eps}
\end{equation}
for any primitive Dirichlet character $\chi$ modulo $q>1$;
see \S\ref{sec:pfThm1}. Since $\EuScript{B}\in C_c^\infty(\R^+)$ is arbitrary,
${\rm GRH}[\frac{9}{10}]$ follows for $L(s,\chi)$.
Conversely, under ${\rm GRH}[\frac{9}{10}]$, we see that \eqref{eq:thorne} holds for every
primitive character $\chi\bmod q$. This leads to similar bounds on
the sums \eqref{eq:einstein} and \eqref{eq:feynman}; see \S\ref{sec:pfThm2}.

\bigskip

{\large\section{Integral bounds}
\label{sec:integralbounds}}

In this section, we establish certain integral bounds that are needed
in the proof of Theorem~\ref{thm:vonMangoldt-twist} (see
\S\ref{sec:vonMangoldt-twist}).

\bigskip

\begin{lemma}\label{lem:mandalorian}
Let $\xi>0$, $\EuScript{B}\in C_c^\infty(\R^+)$, and $X>10$. Then
$$
\int_{\R^+}\Big(1-\frac{1}{u^3-u}\Big)\e(-u\xi)\EuScript{B}(u/X)\,du
\,\mathop{\ll}_{\xi,\EuScript{B}}\,X^{-1},
$$
where the implied constant depends only on $\xi$ and $\EuScript{B}$.
\end{lemma}

\begin{proof}
Put
$$
g(u)\defeq\Big(1-\frac{1}{u^3-u}\Big)\EuScript{B}\Big(\frac{u}{X}\Big)
\qquad(u\in\R^+).
$$
Then $g$ and all of its derivatives are smooth
and supported on the set $X\,\text{\rm supp}(\EuScript{B})$.
In particular, $g''(u)=0$ unless $u\,\mathop{\asymp}_\EuScript{B}\, X$, in
which case one has $g''(u)\,\mathop{\ll}_\EuScript{B}\,X^{-2}$.
Integrating by parts twice, we have (since $g$ and $g'$ vanish at $0$ and $\infty$)
$$
\int_{\R^+}\e(-u\xi)g(u)\,du
=\int_{\R^+}
\frac{\e(-u\xi)}{(-2\pi i\xi)^2}\,g''(u)\,du
\,\mathop{\ll}_\EuScript{B}\,\xi^{-2}X^{-2}
\,{\rm meas}\big(X\,\text{\rm supp}(\EuScript{B})\big),
$$
and the lemma follows.
\end{proof}

The next two technical lemmas are needed in the proof of
Lemma~\ref{lem:the-child} below.

\begin{lemma}\label{lem:technical}
Fix $u_0\in\R^+$, and let $\cL,\cD:C_c^\infty(\R^+)\to C_c^\infty(\R^+)$
be the linear operators defined by
$$
\cL F(u)\defeq\frac{uF(u)}{(u_0-u)}
\mand
\cD F(u)\defeq F'(u)\qquad(u\in\R^+).
$$
For any integer $k\ge 0$, let $(\cD\cL)^k$ and $\cL(\cD\cL)^k$ be the linear
operators given by
$$
[\cD\cL]^k\defeq
\mathop{\underbrace{~\cD\circ\cL\circ\cdots\circ\cD\circ\cL~}}
\limits_{\text{\rm $k$ copies each of $\cD$ and $\cL$}}
\mand
\cL[\cD\cL]^k\defeq\cL\circ[\cD\cL]^k.
$$
Then for every $F\in C_c^\infty(\R^+)$ and every integer $k\ge 0$ we have
\dalign{
[\cD\cL]^kF(u)&\,\mathop{\ll}_{F,k}\,
\max\big\{|u-u_0|^{-k},|u-u_0|^{-2k}\big\},\\
\cL[\cD\cL]^kF(u)&\,\mathop{\ll}_{F,k}\,
\max\big\{|u-u_0|^{-k},|u-u_0|^{-2k-1}\big\},
}
for all $u\in\R^+$, $u\ne u_0$, where the implied constants depend
only on $F$ and $k$.
\end{lemma}

\begin{proof}
Let $F\in C_c^\infty(\R^+)$ be fixed in what follows.
For any integers $A,B\ge 0$ and a real number $C>0$,
let $\EuScript{V}(A,B,C)$ denote the set of functions $G\in C_c^\infty(\R^+)$
of the form
$$
G(u)=\ssum{h,i,j\ge 0\\h\le A\le j\\i+j=h+B}
c_{h,i,j}\frac{u^hF^{(i)}(u)}{(u_0-u)^j},
$$
where the sum runs over nonnegative integers $h,i,j$,
and the coefficients $c_{h,i,j}$ are complex numbers
satisfying $|c_{h,i,j}|\le C$. If $F$ is supported on the
interval $[a,b]$, where $0<a<b<\infty$, then the trivial bound
\begin{equation}\label{eq:eureka0}
\big|G(u)\big|\le C'_{A,B}\,C\max\{b^A,1\}
\max\limits_{0\le i\le A+B}\big|F^{(i)}(u)\big|
\max_{A\le j\le A+B}|u-u_0|^{-j}
\end{equation}
holds for all $G\in\EuScript{V}(A,B,C)$, where
$$
C'_{A,B}\defeq\big|\{h,i,j\ge 0:h\le A\le j,~i+j=h+B\}\big|
\le\frac12(A+1)(A+2B+2).
$$

Next, noting that
$$
\cL G(u)=\ssum{h,i,j\ge 0\\h\le A\le j\\i+j=h+B}c_{h,i,j}
\frac{u^{h+1}F^{(i)}(u)}{(u_0-u)^{j+1}},
$$
it follows that
\begin{equation}\label{eq:eureka1}
\cL:\EuScript{V}(A,B,C)\to\EuScript{V}(A+1,B,C).
\end{equation}
Similarly, since
$$
\cD G(u)=\ssum{h,i,j\ge 0\\h\le A\le j\\i+j=h+B}c_{h,i,j}
\bigg\{\frac{hu^{h-1}F^{(i)}(u)}{(u_0-u)^j}
+\frac{u^hF^{(i+1)}(u)}{(u_0-u)^j}
-\frac{ju^hF^{(i)}(u)}{(u_0-u)^{j+1}}\bigg\},
$$
we see that
\begin{equation}\label{eq:eureka2}
\cD:\EuScript{V}(A,B,C)\to\EuScript{V}\big(A,B+1,(2A+B+1)C\big).
\end{equation}

Observe that $F$ itself lies in $\EuScript{V}(0,0,1)$. Therefore,
an inductive argument using \eqref{eq:eureka1} and \eqref{eq:eureka2} shows that
for every integer $k\ge 0$ we have
$$
[\cD\cL]^kF\in\EuScript{V}(k,k,k!\cdot 3^k)
\mand
\cL[\cD\cL]^kF\in\EuScript{V}(k+1,k,k!\cdot 3^k).
$$
Using \eqref{eq:eureka0}, the result follows.
\end{proof}

\begin{lemma}\label{lem:funnyexpint}
For any positive number $\lambda$, we have
$$
\mathop{\int_\R}_{(|u|>\lambda)}\mathrm{e}^{-iu^2}\,du
=\frac{\mathrm{e}^{-i\lambda^2}}{i\lambda}+i
\mathop{\int_\R}_{(|u|>\lambda)}\mathrm{e}^{-iu^2} u^{-2}\,du
\ll\lambda^{-1}.
$$
\end{lemma}

The next result is our primary tool; its proof is based on the well known
stationary phase method.

\begin{lemma}\label{lem:the-child}
Let $\xi>0$, $\EuScript{B}\in C_c^\infty(\R^+)$, and $X>10$. Suppose
 $\EuScript{B}$ is supported on the interval $[a,b]$, where $0<a<b<\infty$.
For any real number $\gamma\ne 0$, consider the integral $\cI(\gamma)$
defined by
$$
\cI(\gamma)\defeq\int_{\R^+}\e(-u\xi)\EuScript{B}(u/X)u^{-1/2+i\gamma}\,du.
$$
Put
\begin{equation}\label{eq:ast-bst-defn}
a_\star\defeq\frac{1+2a-\sqrt{1+4a}}{2},\qquad
b_\star\defeq\frac{1+2b+\sqrt{1+4b}}{2},\qquad
\gamma_\star\defeq\frac{\gamma}{2\pi\xi X}.
\end{equation}
Then, for $\gamma_\star\not\in[a_\star,b_\star]$ we have
$$
\cI(\gamma)\,\mathop{\ll}_{\xi,\EuScript{B}}\,X^{1/2}
\max\big\{X^{-2}|\gamma|^{-2},|\gamma|^{-4}\big\},
$$
where the implied constant depends only on $\xi$ and $\EuScript{B}$.
On the other hand, for $\gamma_\star\in[a_\star,b_\star]$ we have
$$
\cI(\gamma)=\xi^{-1/2-i\gamma}
\e\Big(\frac{\gamma}{2\pi}\log\frac{\gamma}{2\pi\mathrm{e}}+\frac78\Big)
\EuScript{B}(\gamma_\star)+O_{\xi,\EuScript{B}}(X^{-1/10}),
$$
where the implied constant depends only on $\xi$ and $\EuScript{B}$.
\end{lemma}

\begin{proof}
Making the change of variables $u\mapsto Xu$, we have
\begin{equation}\label{eq:solar1}
\cI(\gamma)=X^{1/2+i\gamma}\int_{\R^+}\mathrm{e}^{iXf(u)}g(u)\,du
=X^{1/2+i\gamma}\,\cJ\quad\text{(say)},
\end{equation}
where
$$
f(u)\defeq -2\pi\xi u+\frac{\gamma\log u}{X},\qquad
g(u)\defeq \frac{\EuScript{B}(u)}{u^{1/2}}.
$$
Note that $\gamma_\star=\gamma/(2\pi\xi X)$
is the only real number for which $f'(\gamma_\star)=0$.

For any given $\Delta>0$, we write
$\cJ=\cJ_\infty+\cJ_\star$
with
$$
\cJ_\infty\defeq\mathop{\int_{\R^+}}_{(|u-\gamma_\star|>\Delta)}
\mathrm{e}^{iXf(u)}g(u)\,du,\qquad
\cJ_\star\defeq\mathop{\int_{\R^+}}_{(|u-\gamma_\star|\le\Delta)}
\mathrm{e}^{iXf(u)}g(u)\,du.
$$
We study $\cJ_\infty$ first.
Let $\cL$ and $\cD$ be the operators defined in
Lemma~\ref{lem:technical} with $u_0\defeq\gamma_\star$,
and put
$$
g_k\defeq [\cD\cL]^kg
\mand
\tilde g_k\defeq \cL[\cD\cL]^kg\qquad(k\ge 0).
$$
According to Lemma~\ref{lem:technical},
\begin{equation}\label{eq:signals}
g_k(u)\,\mathop{\ll}_{\EuScript{B},k}\,
\max\big\{|u-\gamma_\star|^{-k},|u-\gamma_\star|^{-2k}\big\}
\qquad(u\ne\gamma_\star).
\end{equation}
Taking into account that $f'(u)=2\pi\xi u^{-1}(\gamma_\star-u)$, we have
$$
\cJ_\infty=\frac{1}{2\pi\xi}\mathop{\int_{\R^+}}_{(|u-\gamma_\star|>\Delta)}
\mathrm{e}^{iXf(u)}f'(u)\cdot\tilde g_0(u)\,du
=-\frac{1}{2\pi i\xi X}\mathop{\int_{\R^+}}_{(|u-\gamma_\star|>\Delta)}
\mathrm{e}^{iXf(u)} g_1(u)\,du,
$$
where we have used integration by parts in the last step.
Similarly, by induction on $k$, we see that
\begin{equation}\label{eq:signals2}
\cJ_\infty=\frac{1}{(-2\pi i\xi X)^k}\mathop{\int_{\R^+}}_{(|u-\gamma_\star|>\Delta)}
\mathrm{e}^{iXf(u)} g_k(u)\,du.
\end{equation}

The numbers $a_\star$ and $b_\star$ defined in \eqref{eq:ast-bst-defn} have
the property that if $\gamma_\star\not\in[a_\star,b_\star]$, then
$|u-\gamma_\star|>|\gamma_\star|^{1/2}$ for all $u\in[a,b]$. Hence, choosing
$\Delta\defeq|\gamma_\star|^{1/2}$ in the case that
$\gamma_\star\not\in[a_\star,b_\star]$, it follows that $\cJ_\infty=\cJ$
and $\cJ_\star=0$. Setting $k\defeq 4$ and combining \eqref{eq:signals}
and \eqref{eq:signals2}, we derive the bound
$$
\cJ\,\mathop{\ll}_\EuScript{B}\,(\xi X)^{-4}
\max\big\{|\gamma_\star|^{-2},|\gamma_\star|^{-4}\big\}
\ll\max\big\{(\xi X)^{-2}|\gamma|^{-2},|\gamma|^{-4}\big\}
$$
In view of \eqref{eq:solar1}, we obtain the first statement of the lemma.

From now on, we assume that $\gamma_\star\in[a_\star,b_\star]$.
We further assume that
$X$ is large enough (depending on $\EuScript{B}$) so that
\begin{equation}\label{eq:Delta-constraint}
\Delta\defeq X^{-2/5}<\min\{1,2^{-1/2}a_\star\}.
\end{equation}
Write $\cJ=\cJ_\infty+\cJ_\star$ as before with this $\Delta$.
Setting $k\defeq 5$ (say) and combining \eqref{eq:signals}
and \eqref{eq:signals2}, we derive the bound
\begin{equation}\label{eq:solar2}
\cJ_\infty\,\mathop{\ll}_{\xi,\EuScript{B}}\,X^{-1}.
\end{equation}

Turning to the estimate of $\cJ_\star$, observe that the upper bound
\eqref{eq:Delta-constraint} implies that the interval
$\Omega_\Delta\defeq[\gamma_\star-\Delta,\gamma_\star+\Delta]$ lies
entirely inside $\R^+$; in particular,
$$
\cJ_\star\defeq\int_{\Omega_\Delta}\mathrm{e}^{iXf(u)}g(u)\,du.
$$
Since $\gamma_\star\in[a_\star,b_\star]$ (and thus,
$\gamma_\star\,\mathop{\asymp}_\EuScript{B}\,1)$, the estimate
$$
f'''(u)=\frac{4\pi\xi\gamma_\star}{u^3}\,\mathop{\asymp}_\EuScript{B}\,\xi
$$
holds uniformly
for all $u\in\Omega_\Delta$, hence by Taylor's approximation we have
$$
\mathrm{e}^{iXf(u)}=\mathrm{e}^{iX\{f(\gamma_\star)
+\frac12f''(\gamma_\star)(u-\gamma_\star)^2\}}(1+r_1(u))
$$
for some complex function $r_1$ such that
$$
r_1(u)\,\mathop{\ll}_{\xi,\EuScript{B}}\,X|u-\gamma_\star|^3
\qquad(u\in\Omega_\Delta).
$$
We can also write
$$
g(u)=g(\gamma_\star)+r_2(u),
$$
where $r_2$ satisfies the bound
$$
r_2(u)\,\mathop{\ll}_{\xi,\EuScript{B}}\,|u-\gamma_\star|\qquad(u\in\Omega_\Delta).
$$
Therefore, since
\dalign{
\int_{\Omega_\Delta}\big|r_1(u)\big|\,du
&\,\mathop{\ll}_{\xi,\EuScript{B}}\,X\int_{\gamma_\star-\Delta}^{\gamma_\star+\Delta}
|u-\gamma_\star|^3\,du\ll X\Delta^4,\\
\int_{\Omega_\Delta}\big|r_2(u)\big|\,du
&\,\mathop{\ll}_{\xi,\EuScript{B}}\,\int_{\gamma_\star-\Delta}^{\gamma_\star+\Delta}
|u-\gamma_\star|\,du\ll\Delta^2,
}
and $\Delta^2\ll X\Delta^4=X^{-3/5}$ by \eqref{eq:Delta-constraint},
we derive the estimate
\begin{equation}
\label{eq:solar3}
\begin{split}
\cJ_\star
&=\int_{\Omega_\Delta}\mathrm{e}^{iX\{f(\gamma_\star)
+\frac12f''(\gamma_\star)(u-\gamma_\star)^2\}}(1+r_1(u))
\,(g(\gamma_\star)+r_2(u))\,du\\
&=\mathrm{e}^{-i\gamma}\gamma_\star^{-1/2+i\gamma}\EuScript{B}(\gamma_\star)
\int_{\Omega_\Delta}
\mathrm{e}^{-\pi i\xi X\gamma_\star^{-1}(u-\gamma_\star)^2}\,du
+O_{\xi,\EuScript{B}}(X^{-3/5}),
\end{split}
\end{equation}
where in the last step we used the identities
$$
Xf(\gamma_\star)=-\gamma+\gamma\log\gamma_\star,\qquad
f''(\gamma_\star)=-2\pi\xi\gamma_\star^{-1},\qquad
g(\gamma_\star)\defeq \gamma_\star^{-1/2}\EuScript{B}(\gamma_\star).
$$

Next, we extend the range of integration
in the preceding integral to all of $\R$ (with
an acceptable error). Consider the integral
$$
\cK\defeq\mathop{\int_\R}_{(u\not\in\Omega_\Delta)}
\mathrm{e}^{-\pi i\xi X\gamma_\star^{-1}(u-\gamma_\star)^2}\,du.
$$
Making the change of variables $u\mapsto c\,u+\gamma_\star$,
where $c\defeq\sqrt{\gamma_\star/(\pi\xi X)}$, and applying
Lemma~\ref{lem:funnyexpint}, we have
$$
\cK=c\mathop{\int_\R}_{(|u|>\Delta/c)}
\mathrm{e}^{-iu^2}\,du\ll\,c^2\Delta^{-1}\,\mathop{\ll}_{\xi,\EuScript{B}}\,X^{-3/5}.
$$
Using this bound together with \eqref{eq:solar3}, it follows that
$$
\cJ_\star=\mathrm{e}^{-i\gamma}\gamma_\star^{-1/2+i\gamma}\EuScript{B}(\gamma_\star)
\int_\R\mathrm{e}^{-\pi i\xi X\gamma_\star^{-1}(u-\gamma_\star)^2}\,du
+O_{\xi,\EuScript{B}}(X^{-3/5}).
$$

Combining the previous bound with \eqref{eq:solar2}, we have
\begin{equation}\label{eq:solar5}
\cJ=\mathrm{e}^{-i\gamma}\gamma_\star^{-1/2+i\gamma}\EuScript{B}(\gamma_\star)
\int_\R\mathrm{e}^{-\pi i\xi X\gamma_\star^{-1}(u-\gamma_\star)^2}\,du
+O_{\xi,\EuScript{B}}(X^{-3/5}).
\end{equation}
The integral here can be explicitly evaluated:
$$
\int_\R\mathrm{e}^{-\pi i\xi X\gamma_\star^{-1}(u-\gamma_\star)^2}\,du
=\mathrm{e}^{-i\pi/4}\sqrt{\frac{\gamma_\star}{\xi X}}.
$$
Inserting this result into \eqref{eq:solar5}
and recalling \eqref{eq:solar1}, after some
simplification we find that
$$
\cI(\gamma)=\xi^{-1/2-i\gamma}
\e\Big(\frac{\gamma}{2\pi}\log\frac{\gamma}{2\pi\mathrm{e}}-\frac18\Big)
\EuScript{B}(\gamma_\star)
+O_{\xi,\EuScript{B}}(X^{-1/10}),
$$
and the proof is complete.
\end{proof}

{\large\section{Twisting the von Mangoldt function}
\label{sec:vonMangoldt-twist}}

\begin{theorem}\label{thm:vonMangoldt-twist}
Assume RH.
Let $\xi>0$, $\EuScript{B}\in C_c^\infty(\R^+)$, and $X>10$. Then
$$
\sum_{n\ge 1}\Lambda(n)\,\e(-n\xi)\,\EuScript{B}\Big(\frac{n}{X}\Big)
=-\ssum{\rho=\frac12+i\gamma}\xi^{-1/2-i\gamma}
\EuScript{Z}(\rho)\EuScript{B}\Big(\frac{\gamma}{2\pi\xi X}\Big)
+O_{\xi,\EuScript{B}}(X^{9/10}).
$$
\end{theorem}

\begin{proof}
Our goal is to estimate
$$
\sum_{n\ge 1}\Lambda(n)\,\e(-n\xi)\,\EuScript{B}\Big(\frac{n}{X}\Big)
=\sum_{n\ge 1}\Lambda(n)\varphi(n),
$$
where $\varphi(u)\defeq\e(-u\xi)\,\EuScript{B}(u/X)$.
By the explicit formula
(see, e.g., Iwaniec and Kowalski~\cite[Exercise~5, p.~109]{IwaniecKowalski})
we have
\begin{equation}\label{eq:explform}
\sum_{n\ge 1}\Lambda(n)\varphi(n)
=\int_{\R^+}\Big(1-\frac{1}{(u-1)u(u+1)}\Big)\varphi(u)\,du
-\sum_\rho\hat\varphi(\rho),
\end{equation}
where $\hat\varphi$ is the Mellin transform of $\varphi$ given by
$$
\hat\varphi(s)\defeq\int_{\R^+} \varphi(u)u^{s-1}\,du.
$$
Using Lemma~\ref{lem:mandalorian} to bound the integral in \eqref{eq:explform},
we get that
\begin{equation}\label{eq:FcX}
\sum_{n\ge 1}\Lambda(n)\varphi(n)
=-\sum_\rho\hat\varphi(\rho)+O_{\xi,\EuScript{B}}(X^{-1}).
\end{equation}

Next, for any complex zero $\rho=\tfrac12+i\gamma$ of $\zeta(s)$ we have
$$
\hat\varphi(\rho)\defeq\int_{\R^+}
\e(-u\xi)\,\EuScript{B}(u/X)u^{-1/2+i\gamma}\,du,
$$
which is the integral $\cI(\gamma)$ considered in
Lemma~\ref{lem:the-child}. Defining $a,a_\star,a,b_\star,\gamma_\star$
as in Lemma~\ref{lem:the-child}, it follows that
\begin{equation}\label{eq:sydney}
\sum_{\rho=\frac12+i\gamma}\hat\varphi(\rho)
=\Sigma_1+O_{\xi,\EuScript{B}}(X^{9/10}),
\end{equation}
where
$$
\Sigma_1\defeq\ssum{\rho=\frac12+i\gamma\\
\gamma_\star\in[a_\star,b_\star]}\xi^{-1/2-i\gamma}
\e\Big(\frac{\gamma}{2\pi}\log\frac{\gamma}{2\pi\mathrm{e}}+\frac78\Big)
\EuScript{B}(\gamma_\star).
$$
Note that the error term in \eqref{eq:sydney} is a consequence of
the following bounds on the sums that arise
naturally in our application of Lemma~\ref{lem:the-child}:
$$
\ssum{\rho=\frac12+i\gamma\\
\gamma_\star\not\in[a_\star,b_\star]}X^{1/2}
\max\big\{X^{-2}|\gamma|^{-2},|\gamma|^{-4}\big\}
\,\mathop{\ll}_{\xi,\EuScript{B}}\,X^{1/2},
\qquad
\ssum{\rho=\frac12+i\gamma\\
\gamma_\star\in[a_\star,b_\star]}X^{-1/10}
\,\mathop{\ll}_{\xi,\EuScript{B}}\,X^{9/10}.
$$
The condition ``$\gamma_\star\in[a_\star,b_\star]$'' in the above
definition of $\Sigma_1$ is redundant (indeed,
we have $a_\star<a<b<b_\star$ by \eqref{eq:ast-bst-defn}, hence
$\EuScript{B}(\gamma_\star)=0$ if $\gamma_\star\not\in[a_\star,b_\star]$),
and so we can simply write
$$
\Sigma_1=\ssum{\rho=\frac12+i\gamma}\xi^{-1/2-i\gamma}
\e\Big(\frac{\gamma}{2\pi}\log\frac{\gamma}{2\pi\mathrm{e}}+\frac78\Big)
\EuScript{B}\Big(\frac{\gamma}{2\pi\xi X}\Big).
$$
Next, observe that \eqref{eq:NTSTrelation} implies
$$
\e\Big(\frac{T}{2\pi}\log\frac{T}{2\pi\mathrm{e}}+\frac78\Big)
=\e\big(N(T)-S(T)+O(T^{-1})\big)
=\overline{\e\big(S(T)\big)}+O(T^{-1})
$$
provided that $T>0$ is not the ordinate of a zero of $\zeta(s)$
(since $N(T)\in\Z$). Taking the limit as $T\to\gamma$, we get that
$$
\e\Big(\frac{\gamma}{2\pi}\log\frac{\gamma}{2\pi\mathrm{e}}+\frac78\Big)
=\EuScript{Z}(\rho)+O(\gamma^{-1}),
$$
where (as in \S\ref{sec:intro})
$$
\EuScript{Z}(\rho)\defeq\lim\limits_{T\to\gamma}\overline{\e\big(S(T)\big)}.
$$
Thus, up to an acceptable error, we can replace $\Sigma_1$ in
\eqref{eq:sydney} with the quantity
$$
\Sigma_1\defeq\ssum{\rho=\frac12+i\gamma}\xi^{-1/2-i\gamma}
\EuScript{Z}(\rho)\EuScript{B}\Big(\frac{\gamma}{2\pi\xi X}\Big).
$$
The theorem now follows by combining \eqref{eq:FcX}
and \eqref{eq:sydney}.
\end{proof}

{\large\section{Proof of Theorem~\ref{thm:RHvsGRH}}
\label{sec:pfThm1}}

We are now in a position to prove Theorem~\ref{thm:RHvsGRH}.
Observe that \eqref{eq:superbound} and
Theorem~\ref{thm:vonMangoldt-twist} together yield the bound
(under RH)
\begin{equation}\label{eq:wind}
\sum_{n\ge 1}\Lambda(n)\,\e(-n\xi)\,\EuScript{B}(n/X)
\,\mathop{\ll}_{\xi,\EuScript{B},\eps}\,X^{9/10+\eps}
\end{equation}
for every $\xi\in\Q\cap(0,1)$, $\EuScript{B}\in C_c^\infty(\R^+)$, and $\eps>0$, where
the implied constant depends only on
$\xi$, $\EuScript{B}$ and $\eps$.

Since RH is equivalent to GRH in the case of principal characters,
we can assume $\chi$ is nonprincipal. Moreover, if
$\chi$ is induced from a primitive character~$\chi_\star$, then
$L(s,\chi)$ and $L(s,\chi_\star)$ have the same zeros
in the critical strip since
$$
L(s,\chi)=L(s,\chi_*)
\sprod{p\,\mid\,q\\p\,\nmid\,q_*}(1-\chi_*(p)p^{-s}).
$$
Hence, to prove Theorem~\ref{thm:RHvsGRH} it suffices to show
that ${\rm GRH}[\frac{9}{10}]$ holds for any Dirichlet $L$-function
attached to a \emph{primitive} Dirichlet character $\chi$ of modulus $q>1$.
In this situation, the following identity is well known:
\begin{equation}\label{eq:chi-basic}
\chi(n)=\frac{\tau(\chi)}{q}
\sum_{m\bmod q}\overline\chi(m)\e(-mn/q)
\qquad(n\in\Z),
\end{equation}
where the sum runs over any complete set of residue
classes $m\bmod q$, and
$\tau(\chi)$ is the Gauss sum given by
$$
\tau(\chi)\defeq\sum_{n\bmod q}\chi(n)\e(n/q);
$$
see, e.g., Bump~\cite[Chapter 1]{Bump}. In particular, one has
$$
\chi(n)=\frac{\tau(\chi)}{q}
\sum_{m=1}^{q-1}\overline\chi(m)\e(-mn/q)\qquad(n\in\Z),
$$
and therefore
$$
\sum_{n\ge 1}\Lambda(n)\chi(n)\EuScript{B}(n/X)
=\frac{\tau(\chi)}{q}\sum_{m=1}^{q-1}\overline\chi(m)
\sum_{n\ge 1}\Lambda(n)\e(-mn/q)\EuScript{B}(n/X).
$$
Taking $\xi\defeq m/q$ in \eqref{eq:wind}, each inner sum
is $O_{q,\EuScript{B},\eps}(X^{9/10+\eps})$,
hence we conclude that the bound
$$
\sum_{n\ge 1}\Lambda(n)\chi(n)\EuScript{B}(n/X)
\,\mathop{\ll}_{q,\EuScript{B},\eps}\,X^{9/10+\eps}
$$
holds for any $\eps>0$. Since $\EuScript{B}\in C_c^\infty(\R^+)$ is arbitrary,
such bounds are sufficient to conclude that
${\rm GRH}[\frac{9}{10}]$ holds for $L(s,\chi)$.

\bigskip

{\large\section{Proof of Theorem~\ref{thm:RHvsGRH2}}
\label{sec:pfThm2}}

In view of Theorem~\ref{thm:vonMangoldt-twist} it suffices to show
(assuming RH and ${\rm GRH}[\frac{9}{10}]$) that
\begin{equation}\label{eq:moffs}
\sum_{n\ge 1}\Lambda(n)\,\e(-n\xi)\,\EuScript{B}(n/X)
\,\mathop{\ll}_{\xi,\EuScript{B},\eps}\,X^{9/10+\eps}
\end{equation}
holds for all $\xi\in\Q\cap(0,1)$, $\EuScript{B}\in C_c^\infty(\R^+)$, and $\eps>0$, where
the implied constant depends only on
$\xi$, $\EuScript{B}$ and $\eps$. By partial summation, \eqref{eq:moffs}
is a consequence of the bound
\begin{equation}\label{eq:tarkin}
\sum_{n\le X}\Lambda(n)\,\e(-n\xi)
\,\mathop{\ll}_{\xi,\eps}\,X^{9/10+\eps}
\end{equation}
for every $\xi\in\Q\cap(0,1)$ and $\eps>0$, where
the implied constant depends only on
$\xi$ and $\eps$. Fixing $\xi\defeq m/q\in\Q\cap(0,1)$ with $q>1$
and $(m,q)=1$, \eqref{eq:tarkin}  follows from the bound
\begin{equation}\label{eq:gideon}
\sum_{q<p\le X}(\log p)\,\e(-mp/q)
\,\mathop{\ll}_{q,\eps}\,X^{9/10+\eps},
\end{equation}
where the sum runs over primes $p\in(q,X]$. Using Dirichlet characters
to detect a reduced residue class modulo $q$, we have
\dalign{
\sum_{q<p\le X}(\log p)\,\e(-mp/q)
&=\ssum{a\in(\Z/q\Z)^\times}\e(-ma/q)
\ssum{q<p\le X\\p\equiv a\bmod q}\log p\\
&=\frac{1}{\phi(q)}\ssum{a\in(\Z/q\Z)^\times}\e(-ma/q)
\sum_{\chi\bmod q}\overline\chi(a)
\ssum{q<p\le X}(\log p)\,\chi(p),
}
where $\phi$ is the Euler function.
When $(a,q)>1$, we have $\chi(a)=0$ for all characters $\chi\bmod q$; therefore,
we can extend the above sum over $a$ to include all elements of $\Z/q\Z$. Then,
noting that $\ssum{a\in\Z/q\Z}\e(-ma/q)=0$ holds when $m$ and $q$ are coprime,
we see that the principal character $\chi_0\bmod q$ makes no contribution; that is,
$$
\sum_{q<p\le X}(\log p)\,\e(-mp/q)
=\frac{1}{\phi(q)}\ssum{a\in\Z/q\Z}\e(-ma/q)
\ssum{\chi\bmod q\\\chi\ne\chi_0}\overline\chi(a)
\ssum{q<p\le X}(\log p)\,\chi(p).
$$
Finally, for any nonprincipal character $\chi\bmod q$ we have
$$
\ssum{q<p\le X}(\log p)\,\chi(p)
=\ssum{n\le X}\Lambda(n)\chi(n)
+O_{q,\eps}(X^{1/2+\eps})\,\mathop{\ll}_{q,\eps}\,X^{9/10+\eps},
$$
where we have used ${\rm GRH}[\frac{9}{10}]$
for $L(s,\chi)$ in the last step.
Applying this bound in the previous identity and summing over $a$ and $\chi$,
we obtain \eqref{eq:gideon}, which  completes the proof of Theorem~\ref{thm:RHvsGRH2}.

\bigskip

{\large\section{Acknowledgements}}

The author thanks Henryk Iwaniec for some stimulating conversations.


\section{Supplemental Material}

\subsection{Different expressions for the quarkonium correlator present in the literature}

In this section we review the different definitions present in the literature~\cite{Brambilla:2016wgg,Brambilla:2017zei,Yao:2020eqy} that appear in the quantum transport equations for quarkonium. The following discussion will further illuminate the nature of the difference between the correlators~\eqref{HQ-corr} and~\eqref{QA-corr}.

We start from the definition of the chromoelectric correlator for the quarkonium transport equation as given in the main text~\eqref{QA-corr}, which appears in the formulation of the open quantum system for quarkonium in the quantum optical limit~\cite{Yao:2020eqy}:
\begin{eqnarray}
\label{app:QA}
g_E^{\rm Q\bar{Q}}(t) = g^2 T_F \left\langle  \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab}
   {F}_{0i}^b(0) \right) \right\rangle \, .
\end{eqnarray}
In the quantum optical limit, it is the correlator~\eqref{app:QA} at finite frequency that contributes to the quarkonium dissociation and recombination rates:
\begin{eqnarray}
g_E^{Q\bar{Q}}(p_0) = g^2 T_F \int_{-\infty}^{+\infty} \!\! \diff t \, e^{ip_0 t} \left\langle  \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) \right) \right\rangle \,.
\end{eqnarray}
On the other hand, in the quantum Brownian motion limit studied in Refs.~\cite{Brambilla:2016wgg,Brambilla:2017zei}, it is the zero frequency of $g_E^{Q\bar{Q}}$ that matters in the quarkonium transport (see also Ref.~\cite{Yao:2021lus}). However, the expression given in Refs.~\cite{Brambilla:2016wgg,Brambilla:2017zei} (see also Eq.~(2.13) of Ref.~\cite{Brambilla:2020qwo}) is
\begin{align}
\kappa^{Q\bar{Q}} \equiv \frac{g^2T_F}{N_c} {\rm Re} \int_{-\infty}^{+\infty}\!\!\diff t \, \Big\langle \mathcal{T}\Big( \widetilde{E}_i^a(t) \widetilde{E}_i^a(0) \Big) \Big\rangle\,,
\end{align}
where $\mathcal{T}$ denotes time-ordering and $\widetilde{E}_i^a(t) = U_{[-\infty,t]} E_i^a(t) U_{[t,-\infty]} = U_{[-\infty,t]} F_{0i}^a(t) U_{[t,-\infty]}$ with $U$ representing a fundamental Wilson line. The expression of $\kappa^{Q\bar{Q}}$ looks different from our expression $g_E^{Q\bar{Q}}(p_0=0)$ here. In the following, we will show they are equivalent for $p_0=0$. First we find
\begin{align}
\int_{-\infty}^{+\infty} \!\! \diff t \left\langle \mathcal{T} \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) \right) \right\rangle
& = \int_{0}^{+\infty} \!\! \diff t \left\langle \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) \right) \right\rangle + \int_{-\infty}^{0} \!\! \diff t \left\langle \left( {F}_{0i}^a(0) \mathcal{W}_{[0, t]}^{ab}  {F}_{0i}^b(t) \right) \right\rangle \\
& = \int_{0}^{+\infty} \!\! \diff t \left\langle \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) \right) \right\rangle + \int_{0}^{+\infty} \!\! \diff t \left\langle \left( {F}_{0i}^a(0) \mathcal{W}_{[0, -t]}^{ab}  {F}_{0i}^b(-t) \right) \right\rangle \nonumber\\
& = 2 \int_{0}^{+\infty} \!\! \diff t \left\langle \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) \right) \right\rangle \,,\nonumber
\end{align}
where we have relabeled the color indexes $a$ and $b$ in the second term on the first line, flipped the sign of $t$ in the second term on the second line and used translational invariance in $t$ on the last line. Then we can show
\begin{align}
{\rm Re} \int_{-\infty}^{+\infty} \!\! \diff t \left\langle \mathcal{T} \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) \right) \right\rangle & = 2\,{\rm Re} \int_{0}^{+\infty} \!\! \diff t \left\langle \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) \right) \right\rangle \\
& = \int_{0}^{+\infty} \!\! \diff t \left\langle \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) \right) \right\rangle + \int_{0}^{+\infty} \!\! \diff t \left\langle \left( {F}_{0i}^a(0) \mathcal{W}_{[0, t]}^{ab} {F}_{0i}^b(t) \right) \right\rangle \nonumber\\
& = \int_{0}^{+\infty} \!\! \diff t \left\langle \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) \right) \right\rangle + \int_{-\infty}^{0} \!\! \diff t \left\langle \left( {F}_{0i}^a(0) \mathcal{W}_{[0, -t]}^{ab} {F}_{0i}^b(-t) \right) \right\rangle \nonumber\\
& = \int_{0}^{+\infty} \!\! \diff t \left\langle \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) \right) \right\rangle + \int_{-\infty}^{0} \!\! \diff t \left\langle \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) \right) \right\rangle \nonumber\\[4pt]
& \propto g_E^{Q\bar{Q}}(p_0=0) \,, \nonumber
\end{align}
where we have relabeled the color indexes $a$ and $b$ in the second term on the second line, flipped the sign of $t$ in the second term on the third line and used translational invariance in $t$ in the second term on the second-to-last line. Finally we study the relation between $\langle \mathcal{T} ( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) ) \rangle$ and $\langle \mathcal{T} ( \widetilde{E}_i^a(t) \widetilde{E}_i^a(0) ) \rangle$, both of which can be studied using the closed-time path integral methods~\cite{keldysh1965diagram}. Since both correlators are time-ordered, we can insert all the fields contained in the correlators on the time-ordered branch of the Schwinger-Keldysh contour, i.e., the upper branch with type-$1$ fields. Then we can use the standard SU$(N_c)$ Wilson line algebra to derive
\begin{align}
\label{eq:difference-correlators-explicit}
& g^2 T_F \left\langle \mathcal{T} \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab} {F}_{0i}^b(0) \right) \right\rangle \\
&= g^2 T_F \int \! DA_E DA_1 DA_2 \, e^{iS[A_1]-iS[A_2]-S_E[A_E]} F_{0i}^a[A_1](t) \, \mathcal{W}_{[t, 0]}^{ab}[A_1] \, F_{0i}^b[A_1](0) \nonumber \\
    &= g^2 \int \! DA_E DA_1 DA_2 \, e^{iS[A_1]-iS[A_2]-S_E[A_E]} \, {\rm Tr}_{c} \! \left\{ F_{0i}[A_1](t) \, U_{[t, 0]}[A_1] \, F_{0i}[A_1](0) \,  U_{[0, t]}[A_1] \right\}  \nonumber \\
    &= g^2 \int \! DA_E DA_1 DA_2 \, e^{iS[A_1]-iS[A_2]-S_E[A_E]} \, {\rm Tr}_{c} \! \left\{ U_{[-\infty, t]}[A_1] \, F_{0i}[A_1](t) \, U_{[t, 0]}[A_1] \, F_{0i}[A_1](0) \, U_{[0, -\infty]}[A_1] \right\}  \nonumber \\
    &= g^2 \left\langle \mathcal{T} \, {\rm Tr}_c \left( U_{[-\infty,t]} F_{0i}(t) U_{[t,0]} F_{0i}(0) U_{[0,-\infty]} \right) \right\rangle \nonumber \\[4pt]
& = g^2T_F \langle \mathcal{T} ( \widetilde{E}_i^a(t) \widetilde{E}_i^a(0) ) \rangle \,,
 \nonumber
\end{align}
where the subscripts $1$, $2$ and $E$ denote the type-1, type-2 and Euclidean fields.
Putting everything together, we have proved that $\kappa^{Q\bar{Q}}$ defined in Refs.~\cite{Brambilla:2016wgg,Brambilla:2017zei} and our expression $g_E^{Q\bar{Q}}(p_0=0)$ are the same, up to a trivial normalization factor.

We want to emphasize that the second-to-last line of this expression~\eqref{eq:difference-correlators-explicit} does not match the correlator that defines the heavy quark diffusion coefficient $\kappa^Q$~\cite{Casalderrey-Solana:2006fio}, which is given by
\begin{eqnarray}
\kappa^Q \propto g^2 \int_{-\infty}^{+\infty} \!\!\diff t  \left\langle {\rm Tr}_c \left( U_{[-\infty,t]} F_{0i}(t) U_{[t,0]} F_{0i}(0) U_{[0,-\infty]} \right) \right\rangle  = \int_{-\infty}^{+\infty} \!\!\diff t\, g_E^Q(t) \, .
\end{eqnarray}
The key difference between the second-to-last line of Eq.~\eqref{eq:difference-correlators-explicit} and $g_E^Q$ is the operator ordering: in the former case the operators are time-ordered while in the latter they are ordered in the sequence as shown. Conceptually they are different in the sense of Figure~\ref{fig:correlator}: the Wilson loop in~\eqref{HQ-corr} is interrupted by the (thermal) trace over states, whereas~\eqref{QA-corr} can be written in a way such that the Wilson lines only appear between the two chromoelectric field operators. In the original formulation of the heavy quark diffusion coefficient~\cite{Casalderrey-Solana:2006fio}, the Wilson line configuration wraps around the closed-time Schwinger-Keldysh contour with a winding number equal to one. The Euclidean calculation of the heavy quark diffusion coefficient~\cite{Burnier:2010rp} also has this feature (see~\cite{Eller:2019spw} for an explicit proof that the Minkowski formulation~\cite{Casalderrey-Solana:2006fio} and the Euclidean formulation~\cite{Burnier:2010rp} give the same result). The Wilson line configuration in the correlator for quarkonium has a winding number equal to zero.
This mathematical difference has physical origin as discussed in the main text.
Therefore, these two quantities $\kappa^{Q\bar{Q}}$ and $\kappa^Q$ (or more generally $g_E^{\rm Q\bar{Q}}$ and $g_E^Q$) cannot be used interchangeably.

The first verification that $g_E^{\rm Q\bar{Q}}$ and $g_E^Q$ are different was achieved in~\cite{Burnier:2010rp}, which can be seen by comparing the results obtained therein for $g_E^Q$ to the results of~\cite{Eidemuller:1997bb}, which first computed the correlator $g_E^{\rm Q\bar{Q}}$ in vacuum. This difference was subsequently verified in~\cite{Binder:2021otw}, which calculates $g_E^{\rm Q\bar{Q}}$ both in vacuum and at finite temperature.
In the present paper we further study their difference with a more general gauge choice, and show the breakdown of a naive axial gauge calculation.

\subsection{Chromoelectric correlator in axial gauge}
Here we study the chromoelectric field correlators for heavy quarks and quarkonia in axial gauge, in which the Wilson lines become identities and thus can be neglected. The time-ordered chromoelectric field correlator is
\begin{eqnarray}
g_{E,\,T}^{\rm Axial}(p_0) = g^2 \int_{-\infty}^{+\infty} \! \diff t\, e^{i p_0 t} \langle 0 |  \mathcal{T}\{ {E}^a_{i}(t,{\boldsymbol x})  {E}^a_{i}(0,{\boldsymbol x}) \} | 0 \rangle\,,
\end{eqnarray}
and we want to calculate it in axial gauge at next-to-leading order (NLO). We will focus on the gluon polarization diagram, which contributes to the correlator at NLO. First we work out the gluon propagator in axial gauge. The free part of the gauge boson Lagrangian plus the gauge-fixing term in momentum space can be written as
\begin{align}
\frac{i}{2} \int \diff^4 k \, A^{\mu a}(-k) \Big( - g_{\mu\nu}(k^2 +i\varepsilon) +k_\mu k_\nu - \frac{1}{\xi}n_\mu n_\nu \Big) A^{\nu a}(k) \,,
\end{align}
where $\varepsilon$ comes from the boundary condition of the path integral at $t=\pm\infty$ and $\xi$ is a gauge-fixing parameter to be set later. Inverting $i g_{\mu\nu}(k^2 +i\varepsilon) -i k_\mu k_\nu +i n_\mu n_\nu/\xi$ gives the time-ordered gluon propagator
\begin{align}
[D_T(k)]_{\mu\nu}^{ab} = \frac{i\delta^{ab}}{k^2+i\epsilon} \left[ -g_{\mu\nu} + \frac{ n\cdot k(k_\mu n_\nu + k_\nu n_\mu) - [\xi(k^2+i\varepsilon)+t^2] k_\mu k_\nu +i\varepsilon n_\mu n_\nu }{i\varepsilon[\xi(k^2+i\varepsilon)+n^2] + (n\cdot k)^2} \right] \,.
\end{align}
Setting $\xi=0$ and $n_{\mu}=(1,0,0,0)$ for temporal axial gauge and neglecting terms proportional to $\varepsilon$ in the numerator lead to
\begin{align}
[D_T(k)]_{\mu\nu}^{ab} = \frac{i\delta^{ab} P_{\mu\nu}(k) }{k^2+i\epsilon} \equiv D_T(k) \delta^{ab} P_{\mu\nu}(k) \,,
\end{align}
where
\begin{align}
    P_{\mu \nu}(k) = -g_{\mu \nu} + \frac{k_0(k_\mu n_\nu + n_\mu k_\nu)}{k_0^2 + i\varepsilon} - \frac{k_\mu k_\nu}{k_0^2 + i\varepsilon} \, .
\end{align}
Then, using the Feynman rules of non-Abelian gauge theory we find the contribution of the gluon polarization diagrams (with its two external legs connected with the two chromoelectric fields) to the time-ordered chromoelectric correlator
\begin{align}
\left. g_{E,\,T}^{\rm Axial}(p_0) \right|_{\rm NLO} &= g^2 \int_{\boldsymbol p} D_T(p)^2 (i p_0 g_i^{\,\sigma'} - i p_i g_0^{\,\sigma'}) P_{\sigma' \sigma}(p) (- i p_0 g_i^{\,\rho'} + i p_i g_0^{\,\rho'}) P_{\rho' \rho}(p) \nonumber \\
& \quad \frac12 \bigg( \delta^{cd} \int_k D_T(k) \delta^{ab} P_{\mu \nu}(k) (-i g^2) \left[ f^{abe} f^{cde} (g^{\mu \rho} g^{\nu \sigma} - g^{\mu \sigma} g^{\nu \rho}) + f^{ace} f^{dbe} (g^{\mu \sigma} g^{\nu \rho} - g^{\mu \nu} g^{\rho \sigma}  ) \right.  \nonumber \\
    & \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad + \left. f^{ade} f^{bce} (g^{\mu \nu} g^{\rho \sigma} - g^{\mu \rho} g^{\nu \sigma}) \right] \nonumber \\
& + \delta^{cc'} \int_k D_T(k) D_T(p-k) \delta^{aa'} \delta^{bb'} P_{\mu\mu'}(k) P_{\nu\nu'}(p-k) g f^{abc} \left[ g^{\mu \nu} (p - 2k)^{\sigma} + g^{\nu \sigma} (k-2p)^{\mu} + g^{\sigma \mu}(p+k)^{\nu} \right] \nonumber \\
    & \quad \quad \quad \quad \quad \quad \quad \quad \times g f^{a'b'c'} \left[ g^{\mu' \nu'} (2k-p)^{\rho} + g^{\nu' \rho} (2p-k)^{\mu'} + g^{\rho \mu'}(-p-k)^{\nu'} \right] \bigg) \,,
\end{align}
where $\int_{\boldsymbol p} \equiv \int\frac{\diff^d p}{(2\pi)^d}$ and $\int_k \equiv \int\frac{\diff^D k}{(2\pi)^D}$.
Our strategy to evaluate these integrals is to do the ${\boldsymbol p}$ and ${\boldsymbol k}$ integrals first, using dimensional regularization in $d = 3 - \tilde{\epsilon}$ spatial dimensions for both of them. $D = 4 - \tilde{\epsilon}$ is the total number of spacetime dimensions. (The calculation is only consistent if we use the same dimensionality for both ${\boldsymbol p}$ and ${\boldsymbol k}$ integrals.) We leave the integral over $k_0$ to be done at the end of the calculation.
We proceed by reducing the integral into a handful of integral structures $\tilde{I}_i(p,k)$ and their respective numerators $N_i(p_0,k_0)$ that do not depend on spatial momenta, which gives
\begin{eqnarray}
g_{E,\,T}^{\rm Axial}(p_0) = \frac{N_c (N_c^2-1) g^4 p_0^2 }{2} \int_{-\infty}^{+\infty} \! \frac{\diff k_0} {2\pi} \int \frac{\diff^d {\boldsymbol p} \diff^d {\boldsymbol k} }{(2\pi)^{2d}} \sum_i N_i(p_0,k_0) \tilde{I}_i(p,k) \, .
\end{eqnarray}
Also, we denote $I_i(p_0,k_0) = \int \frac{\diff^d {\boldsymbol p} \diff^d {\boldsymbol k} }{(2\pi)^{2d}} \tilde{I}_i(p,k)$.

Below we list the resulting integral structures, accompanied by their respective numerators:
\begin{enumerate}
    \item \begin{align}
    I_1 &= \int_{\boldsymbol p} \int_{\boldsymbol k} \frac{1}{(p^2+i\varepsilon)^2 (k^2 + i\varepsilon) ((p-k)^2 + i\varepsilon) } \, , \\
    N_1 &= 0 \, , \end{align}
    \item \begin{align}
    I_2 &= \int_{\boldsymbol p} \int_{\boldsymbol k} \frac{1}{(p^2+i\varepsilon) (k^2 + i\varepsilon) ((p-k)^2 + i\varepsilon) } \nonumber \\ &= - \frac{\Gamma(4-D)}{(4\pi)^{D-1} } \int_0^1 \diff x \diff y \frac{ \left[ - y (k_0 - x p_0)^2 - (1 - y + y x(1-x)) p_0^2 - i\varepsilon \right]^{D-4} }{y^{\frac{D-3}{2}} (1-y+yx(1-x))^{\frac{D-1}{2}}} \, , \\
    N_2 &=  \frac{4 (D-2) (k_0^2 - k_0 p_0 + p_0^2)^2 }{k_0 (k_0 - p_0) p_0^2 } \, , \end{align}
    \item \begin{align}
    I_3 &= \int_{\boldsymbol p} \int_{\boldsymbol k} \frac{1}{(p^2+i\varepsilon)^2 (k^2 + i\varepsilon)} = - (-i)^{2D-8} \frac{\Gamma \! \left( \frac{3-D}{2} \right) \Gamma \! \left( \frac{5-D}{2} \right) }{(4\pi)^{D-1}} \frac{|k_0 p_0|^{D-3}}{p_0^2} \, , \\
    N_3 &=  \frac{2 (D-2) \left[ 2(D-1) k_0^2 - (D-2) p_0 k_0 + 2(D-1) p_0^2 \right] }{(D-1) k_0 p_0 } \, , \end{align}
    \item \begin{align}
    I_4 &= \int_{\boldsymbol p} \int_{\boldsymbol k} \frac{1}{(p^2+i\varepsilon) (k^2 + i\varepsilon)}  = (-i)^{2D-6} \frac{\Gamma\!\left(\frac{3-D}{2} \right)^2}{(4\pi)^{D-1}} |k_0 p_0|^{D-3} \, , \\
    N_4 &=  \frac{2 \left[ 2(D-2) k_0^4 - 3(D-2) k_0^3 p_0 + D\, k_0^2 p_0^2 - (D-2) k_0 p_0^3 + (D-2) p_0^4 \right]  }{ k_0^2 p_0^3 (k_0 - p_0) } \, , \end{align}
    \item \begin{align}
    I_5 &= \int_{\boldsymbol p} \int_{\boldsymbol k} \frac{1}{((p-k)^2+i\varepsilon) (k^2 + i\varepsilon)} = (-i)^{2D-6} \frac{\Gamma\!\left(\frac{3-D}{2} \right)^2}{(4\pi)^{D-1}} |k_0 (p_0-k_0) |^{D-3} \, , \\
    N_5 &=  \frac{D-2}{k_0^2} + \frac{D-2}{(k_0-p_0)^2} + \frac{2}{p_0^2} \, , \end{align}
\end{enumerate}
All the other integral structures give vanishing contributions in dimensional regularization. Note that the term 1.~vanishes because the numerator happens to be zero, and the term 3.~also vanishes upon integration over $k_0$ because the integrand is just a polynomial in $k_0$. (In dimensional regularization the limit $\tilde{\epsilon} \to 0$ cannot be taken before performing all integrals that involve $d$, which means the limit should be taken after the $k_0$ integral.) Then, one can show that
\begin{enumerate}
    \item \begin{align}
    \frac{N_c (N_c^2-1) g^4 p_0^2 }{2}  \int_{-\infty}^{+\infty} \frac{\diff k_0}{2\pi}  N_1(p_0,k_0) I_1(p,k) = 0
     \, , \end{align}
    \item \begin{align}
    \frac{N_c (N_c^2-1) g^4 p_0^2 }{2}  \int_{-\infty}^{+\infty} \frac{\diff k_0}{2\pi}  N_2(p_0,k_0) I_2(p,k) = \frac{N_c(N_c^2-1)g^4 p_0^3}{(2\pi)^3} \left[ \frac{11}{6 \tilde{\epsilon} } + \frac{11}{6} \ln \! \left( \frac{\mu^2}{4p_0^2} \right) + \frac{167}{36} + \frac{\pi^2}{3} \right]
     \, , \end{align}
    \item \begin{align}
    \frac{N_c (N_c^2-1) g^4 p_0^2 }{2} \int_{-\infty}^{+\infty} \frac{\diff k_0}{2\pi}  N_3(p_0,k_0) I_3(p,k) = 0
     \, , \end{align}
    \item \begin{align}
    \frac{N_c (N_c^2-1) g^4 p_0^2 }{2} \int_{-\infty}^{+\infty} \frac{\diff k_0}{2\pi}  \left[ N_4(p_0,k_0) I_4(p,k) + N_5(p_0,k_0) I_5(p_0,k_0) \right] = \frac{N_c (N_c^2 - 1) g^4 p_0^3}{ (2\pi)^3} \frac5{12}
     \, . \end{align}
\end{enumerate}

The final contribution to evaluate is from the coupling constant counterterm, since the definition of the chromoelectric correlator contains $g^2$. The contribution reads
\begin{align}
    & (Z_g - 1 ) (N_c^2 - 1) g^2 \int_{\boldsymbol p}  D_T(p) (i p_0 g_i^{\,\sigma} - i p_i g_0^{\,\sigma}) P_{\sigma \sigma'}(p) (- i p_0 g_i^{\,\sigma'} + i p_i g_0^{\,\sigma'}) \nonumber \\
    &= \frac{g^4}{8\pi^2 (D-4)} \frac{11}{3} N_c (N_c^2 - 1) p_0^2 \int_{\boldsymbol p} \frac{i (\delta_{ii} - {\boldsymbol p}^2/p_0^2 )}{p_0^2 - {\boldsymbol p}^2 + i\varepsilon} \nonumber \\
    &= \frac{g^4 N_c (N_c^2 - 1)}{(2\pi)^3} p_0^2 \left[ \frac{11}{24 \pi^2 (D-4)} \frac{1}{(2\pi)^{D-4}} \pi \frac{1}{2p_0} \Omega_{D-1} ( D-2 ) p_0^{D-2}  \right] \nonumber \\
    &= \frac{g^4 N_c (N_c^2 - 1)}{(2\pi)^3} p_0^3 \frac{11}{48\pi} \left[ \frac{1}{ (D-4)} \frac{\Omega_{D-1} ( D-2 ) (p_0/\tilde{\mu} )^{D-4}}{(2\pi)^{D-4}}   \right] \nonumber \\
    &= \frac{g^4 N_c (N_c^2 - 1)}{(2\pi)^3} p_0^3 \left[ - \frac{11}{6\tilde{\epsilon}} + \frac{11}{48\pi} \frac{\partial}{\partial D} \left( \frac{\Omega_{D-1} ( D-2 ) (p_0/\tilde{\mu} )^{D-4}}{(2\pi)^{D-4}} \right)_{D=4} \right] \nonumber \\
    &= \frac{g^4 N_c (N_c^2 - 1)}{(2\pi)^3} p_0^3 \left[ - \frac{11}{6\tilde{\epsilon}} + \frac{11}{12} \left( 1 + \ln( \pi (p_0/(2\pi \tilde{\mu}))^2 ) - ( 2 - \gamma_E - 2\ln(2) ) \right) \right] \nonumber \\
    &=  \frac{g^4 N_c (N_c^2 - 1)}{(2\pi)^3} p_0^3 \left[ - \frac{11}{6\tilde{\epsilon}} - \frac{11}{12} - \frac{11}{12} \ln \! \left( \frac{\mu^2}{4p_0^2} \right)  \right] \,,
\end{align}
where $\mu^2 = 4\pi e^{-\gamma_E} \tilde{\mu}^2$.
Adding everything up, one obtains
\begin{equation}
    \left. g_{E,\,T}^{\rm Axial}(p_0) \right|_{\rm NLO}  = \frac{g^4 N_c (N_c^2 - 1)}{(2\pi)^3} p_0^3 \left[ \frac{11}{12}\ln \! \left( \frac{\mu^2}{4p_0^2} \right) + \frac{149}{36} + \frac{\pi^2}{3} \right],
\end{equation}
which is what we give in the main text.

The corresponding Euclidean correlator can be evaluated in exactly the same way. In the same notation (but with the understanding that $k_0$ and $p_0$ are now Euclidean momenta), the relevant integral structures are
\begin{enumerate}
    \item \begin{align}
    I_1 &= \int_{\boldsymbol p} \int_{\boldsymbol k} \frac{1}{(p^2)^2 k^2 (p-k)^2 } \, , \\
    N_1 &= 0 \, , \end{align}
    \item \begin{align}
    I_2 &= \int_{\boldsymbol p} \int_{\boldsymbol k} \frac{1}{p^2 k^2  (p-k)^2 } \nonumber \\ &=  \frac{\Gamma(4-D)}{(4\pi)^{D-1} } \int_0^1 \diff x \diff y \frac{ \left[  y (k_0 - x p_0)^2 + (1 - y + y x(1-x)) p_0^2 \right]^{D-4} }{y^{\frac{D-3}{2}} (1-y+yx(1-x))^{\frac{D-1}{2}}}  \, , \\
    N_2 &= - \frac{4 (D-2) (k_0^2 - k_0 p_0 + p_0^2)^2 }{k_0 (k_0 - p_0) p_0^2 } \, , \end{align}
    \item \begin{align}
    I_3 &= \int_{\boldsymbol p} \int_{\boldsymbol k} \frac{1}{(p^2)^2 k^2} = \frac{\Gamma \! \left( \frac{3-D}{2} \right) \Gamma \! \left( \frac{5-D}{2} \right) }{(4\pi)^{D-1}} \frac{|k_0 p_0|^{D-3}}{p_0^2} \, , \\
    N_3 &= - \frac{2 (D-2) \left[ 2(D-1) k_0^2 - (D-2) p_0 k_0 + 2(D-1) p_0^2 \right] }{(D-1) k_0 p_0 } \, , \end{align}
    \item \begin{align}
    I_4 &= \int_{\boldsymbol p} \int_{\boldsymbol k} \frac{1}{p^2 k^2 }  =  \frac{\Gamma\!\left(\frac{3-D}{2} \right)^2}{(4\pi)^{D-1}} |k_0 p_0|^{D-3} \, , \\
    N_4 &=  -\frac{2 \left[ 2(D-2) k_0^4 - 3(D-2) k_0^3 p_0 + D k_0^2 p_0^2 - (D-2) k_0 p_0^3 + (D-2) p_0^4 \right]  }{ k_0^2 p_0^3 (k_0 - p_0) } \, , \end{align}
    \item \begin{align}
    I_5 &= \int_{\boldsymbol p} \int_{\boldsymbol k} \frac{1}{(p-k)^2 k^2 } = \frac{\Gamma\!\left(\frac{3-D}{2} \right)^2}{(4\pi)^{D-1}} |k_0 (p_0-k_0) |^{D-3} \, , \\
    N_5 &= - \frac{D-2}{k_0^2} - \frac{D-2}{(k_0-p_0)^2} - \frac{2}{p_0^2} \, , \end{align}
\end{enumerate}
all of which give the same contributions to the correlator as in the time-ordered case.

\subsection{Difference between heavy quark and quarkonium correlators in mixed axial-Feynman gauge}

In this section we show how to calculate the difference between the two types of correlators discussed in the main text in mixed axial-Feynman gauge. To enforce the different operator orderings in the correlation functions, a rigorous calculation of their difference (in particular, the heavy quark correlator) requires introducing the Schwinger-Keldysh contour~\cite{keldysh1965diagram}. As mentioned in the main text, we use a mixed gauge-fixing condition given by
\begin{eqnarray}
G_M^a[A] = \frac{1}{\lambda} n^\mu A_\mu^a (x) + \partial^\mu A^a_\mu(x) \, ,
\end{eqnarray}
over which we perform the standard average over field configurations $\delta(G_M^a[A] - \omega^a(x))$ weighted by a function $\exp(-\frac{i}{2\xi}\int\diff^4x \omega^a\omega^a)$ with a smearing parameter $\xi$~\cite{Srednicki:2007qs,Schwartz:2014sze}. When $\lambda \to \infty$, one recovers Feynman gauge by setting $\xi = 1$ while axial gauge is recovered for $\lambda \to 0$ for any $\xi$. As such, we shall choose $\xi = 1$ throughout. After performing the gauge-fixing procedure in the path integral, the gluon propagators in this gauge are obtained by inverting the kinetic term in the action for the gauge field:
\begin{align}
    &i S_{\rm kin}[A] = \\
    & -\frac12 \int_k \begin{pmatrix} \\ A_{(1)}^{a \nu}(-k) \\ \\ A_{(2)}^{a \nu}(-k) \\ {} \end{pmatrix}^{T} \begin{pmatrix} \begin{matrix} i \Big[ \left(k^2 + i\varepsilon \right) g_{\mu \nu} - \left(1 - \frac{1}{\xi} \right) k_\mu k_\nu \\ \,\,\, + \frac{1}{\xi \lambda^2} n_\mu n_{\nu} + \frac{1}{\xi \lambda} \left( - i k_\mu n_\nu + i k_\nu n_\mu \right)  \Big] \end{matrix}  & 2\varepsilon g_{\mu \nu} \Theta(-k_0) \\ 2\varepsilon g_{\mu \nu} \Theta(-k_0) &  \begin{matrix} -i \Big[ \left(k^2 - i\varepsilon \right) g_{\mu \nu} - \left(1 - \frac{1}{\xi} \right) k_\mu k_\nu \\ \,\,\, + \frac{1}{\xi \lambda^2} n_\mu n_{\nu} + \frac{1}{\xi \lambda} \left( - i k_\mu n_\nu + i k_\nu n_\mu \right)  \Big] \end{matrix} \end{pmatrix} \begin{pmatrix} \\ A_{(1)}^{a \mu}(k) \\ \\ A_{(2)}^{a \mu}(k) \\ {} \end{pmatrix} \nonumber
\end{align}
and we will set $\xi = 1$ in what follows.
The resulting propagators read
\begin{align}
    [D_T(k)]_{\mu\nu}^{ab} &= \frac{i\delta^{ab}}{k^2+i\varepsilon} \left[ -g_{\mu\nu} + \frac{k_\mu n_\nu f + k_\nu n_\mu f^* - n^2 k_\mu k_\nu }{ |f|^2 + i h \varepsilon } \right] \,, \\
    [D_{\bar{T}}(k)]_{\mu \nu}^{ab} &= \frac{-i\delta^{ab}}{k^2-i\varepsilon} \left[ -g_{\mu\nu} + \frac{k_\mu n_\nu f + k_\nu n_\mu f^* - n^2 k_\mu k_\nu }{ |f|^2 - i h \varepsilon } \right] \,, \\
    [D_>(k)]_{\mu\nu}^{ab} &= \frac{2\varepsilon \Theta(k_0) }{(k^2)^2 + \varepsilon^2} \left[ - g_{\mu \nu} + \frac{k_\mu n_\nu F + k_\nu n_\mu F^* -  k_\mu k_\nu H }{|f|^4 + h^2 \varepsilon^2 } \right] \,, \\
    [D_<(k)]_{\mu\nu}^{ab} &= \frac{2\varepsilon \Theta(-k_0) }{(k^2)^2 + \varepsilon^2} \left[ - g_{\mu \nu} + \frac{k_\mu n_\nu F + k_\nu n_\mu F^* -  k_\mu k_\nu H }{|f|^4 + h^2 \varepsilon^2 } \right] \,,
\end{align}
where we have denoted for brevity
\begin{align}
    f &= (n \! \cdot \! k - i\lambda k^2) \,,  \\
    h &= (n^2+ 2\lambda^2 k^2) \,, \\
    F &= k^2 n^2 f + 2 f |f|^2 -  (n \! \cdot \! k) f^2  \,,\\
    H &= 3|f|^2 + (n^2)^2 k^2 - (n \! \cdot \! k) (f + f^*) n^2 \, ,
\end{align}
and ${}^*$ denotes complex conjugation. We have also dropped any terms that have more powers of $\varepsilon$ than listed in the numerators.

We then proceed to evaluate the difference between the correlators
\begin{align}
g_E^{\rm Q}(t) &= g^2 \left\langle {\rm Tr}_c \left( U_{[-\infty,t]} F_{0i}(t) U_{[t,0]} F_{0i}(0) U_{[0,-\infty]} \right) \right\rangle = g^2 \left\langle {\rm Tr}_c \left( U_{[-\infty,t]} F_{0i}(t) U_{[t,+\infty]} U_{[+\infty,0]} F_{0i}(0) U_{[0,-\infty]} \right) \right\rangle \,, \\
g_E^{\rm Q\bar{Q}}(t) &= g^2 T_F \left\langle  \left( {F}_{0i}^a(t) \mathcal{W}_{[t, 0]}^{ab}
   {F}_{0i}^b(0) \right) \right\rangle = g^2 T_F \left\langle  \left( {F}_{0i}^a(t) \mathcal{W}_{[t, +\infty]}^{ac} \mathcal{W}_{[+\infty, 0]}^{cb}
   {F}_{0i}^b(0) \right) \right\rangle \,,
\end{align}
where the last arrangement of operators in each correlator is the one that is better suited for direct evaluation on the Schwinger-Keldysh contour, because in each case we have anti-time ordered operators grouped together on the left and time-ordered operators grouped on the right regardless of whether the time $t$ satisfies $t>0$ or $t<0$. Concretely, in both correlators, we take operators on the left of the Wilson lines extending from $t$ to $t=+\infty$ (including these Wilson lines) to be of type $2$, and operators on the right to be of type $1$.

\begin{figure}
    \begin{subfigure}[t]{0.49\textwidth}
        \centering
        \includegraphics[height=2.2in]{QQbar_EE.pdf}
        \caption{Quarkonium.}
        \label{fig:QQbar_EE}
    \end{subfigure}%
    ~
    \begin{subfigure}[t]{0.49\textwidth}
        \centering
        \includegraphics[height=2.2in]{Q_EE.pdf}
        \caption{Heavy quark.}
        \label{fig:Q_EE}
    \end{subfigure}%
\caption{Feynman diagrams relevant for the difference between the chromoelectric field correlators for quarkonia (left) and heavy quarks (right). The blobs represent the chromoelectric fields while the double/single-dashed lines indicate the adjoint/fundamental Wilson lines. Similar diagrams where the gluon lines originating from the Wilson lines on the right are also included in the calculation.}
\label{fig:5}
\end{figure}

We perform a perturbative expansion on the coupling constant $g$ to calculate the difference
\begin{eqnarray}
g_E^{\rm Q\bar{Q}}(p_0) -  g_E^{\rm Q}(p_0) = \int_{-\infty}^{+\infty} \!\! \diff t \, e^{i p_0 t} \! \left(  g_E^{\rm Q\bar{Q}}(t) -  g_E^{\rm Q}(t) \right) \, .
\end{eqnarray}
In the language of Ref.~\cite{Binder:2021otw}, the difference comes solely from the diagrams of type $(5)$ and $(5r)$, i.e., diagrams with a triple-gauge boson vertex where only one of the three gluon lines is attached to a Wilson line. In the language of Refs.~\cite{Eller:2019spw,Burnier:2010rp}, the difference comes from the diagrams labeled as (j). These diagrams are shown in Fig.~\ref{fig:5}. Here we calculate the correlators directly as defined, rather than the time-ordered version as in the naive axial gauge calculation shown in the previous section. Following the calculation details of Ref.~\cite{Binder:2021otw}, we find these diagrams give
\begin{align}
    g_E^{\rm Q\bar{Q}}(p_0) -  g_E^{\rm Q}(p_0) = \int_{{\boldsymbol p}, k} & T_F g^4 N_c (N_c^2 - 1)  \pi \delta(k_0)  \big[ g_{\mu \nu} (p - 2k)_\rho + g_{\nu \rho} (k - 2p)_{\mu} + g_{\rho \mu} (p+k)_\nu \big] \nonumber \\
    \times &  (p_0 g_{i\rho'} - p_i g_{0\rho'} ) \big((p_0-k_0) g_{i\nu'} - (p_i - k_i) g_{0\nu} \big) \nonumber \\
    \times &  \Big( [D_{>}(p)]^{\rho' \rho} [D_T(p-k)]^{\nu \nu'} [D_T(k)]^{\mu 0} - [D_{\bar{T}}(p)]^{\rho' \rho} [D_>(p-k)]^{\nu \nu'} [D_>(k)]^{\mu 0} \nonumber \\
    & \,\,\,  - [D_{T}(p)]^{\rho \rho' } [D_>(p-k)]^{\nu' \nu } [D_>(k)]^{ 0 \mu}  + [D_{>}(p)]^{\rho \rho' } [D_{\bar{T}}(p-k)]^{\nu' \nu } [D_{\bar{T}}(k)]^{0 \mu} \Big) \, .
\end{align}
Here we have taken the color contractions out of the propagators and therefore only spacetime indices remain. The Dirac delta $\delta(k_0)$ appears when taking the difference between contributions with different pole prescriptions. These prescriptions are determined by the Wilson line regulator $\eta$ (see main text around~\eqref{eq:Wilson-line-regulated}) in each of the two correlators. In fact, this is the only difference between the expressions resulting from the Feynman rules for each diagram. It results in contributions of the form
\begin{eqnarray}
\frac{1}{k_0 + i\eta} - \frac12 \left( \frac{1}{k_0 + i \eta} + \frac{1}{k_0 - i \eta} \right) = \frac12 \left(\frac{1}{k_0 + i \eta} - \frac{1}{k_0 - i \eta} \right) \overset{\eta \to 0}{\longrightarrow} - i \pi \delta(k_0) \nonumber \, ,
\end{eqnarray}
where the first term can be traced back to an adjoint Wilson line attached to $t = +\infty$ (as in the quarkonium correlator), and the second term is the sum of the contributions from a fundamental Wilson line attached to $t=+\infty$ and another one attached to $t=-\infty$ (as in the heavy quark correlator).

The rest of the calculation is tedious, but straightforward. The most sensitive terms come from the propagator that connects the triple-gauge boson vertex with the Wilson line. We list them here explicitly:
\begin{align}
    \left. [D_T(k)]^{\mu 0} \right|_{k_0=0} &= \frac{-i}{{\boldsymbol k}^2-i\varepsilon} \left[ -g^{\mu 0} + \frac{i\lambda {\boldsymbol k}{}^2 k^{\mu} }{\lambda^2 {({\boldsymbol k}{}^2)}{}^2 + (1 - 2\lambda^2 {\boldsymbol k}{}^2 ) i \varepsilon } \right] \,, \\
    \left. [D_{\bar{T}}(k)]^{0 \mu} \right|_{k_0=0} &= \frac{i}{{\boldsymbol k}^2+i\varepsilon} \left[ -g^{ 0 \mu} + \frac{-i\lambda {\boldsymbol k}{}^2 k^{\mu} }{\lambda^2 {({\boldsymbol k}{}^2)}{}^2 - (1 - 2\lambda^2 {\boldsymbol k}{}^2 ) i \varepsilon } \right]  \,, \\
    \left. [D_>(k)]^{\mu 0} \right|_{k_0=0} &= \frac{\varepsilon}{{({\boldsymbol k}{}^2)}{}^2 + \varepsilon^2}  \left[ -g^{\mu 0} + \frac{i \big[ 2\lambda^3 {({\boldsymbol k}{}^2)}{}^3 - \lambda {({\boldsymbol k}{}^2)}{}^2 \big] k^{\mu} }{\lambda^4 {({\boldsymbol k}{}^2)}{}^4 + (1 - 2\lambda^2 {\boldsymbol k}{}^2 ){}^2 \varepsilon^2 } \right]  \,,\\
    \left. [D_>(k)]^{0 \mu} \right|_{k_0=0} &= \frac{\varepsilon}{{({\boldsymbol k}{}^2)}{}^2 + \varepsilon^2}  \left[ -g^{0 \mu} + \frac{-i \big[2\lambda^3 {({\boldsymbol k}{}^2)}{}^3 - \lambda {({\boldsymbol k}{}^2)}{}^2 \big] k^{\mu} }{\lambda^4 {({\boldsymbol k}{}^2)}{}^4 + (1 - 2\lambda^2 {\boldsymbol k}{}^2 ){}^2 \varepsilon^2 } \right] \, .
\end{align}
Performing the index contractions leads to various integral structures. After isolating the contributions where the ${\boldsymbol k}$ momentum flowing in the propagators decouples from ${\boldsymbol p}$, using $\int \frac{\diff^d {\boldsymbol k} }{(2\pi)^d} {(({\boldsymbol k})^{2})}^{n} = 0$  for any integer $n$ in dimensional regularization for $d = 3 -\tilde{\epsilon}$, and using the symmetries of the integrand, one can reduce the expression for the difference to
\begin{align}
    g_E^{\rm Q\bar{Q}}(p_0) -  g_E^{\rm Q}(p_0) = T_F g^4 N_c (N_c^2-1) \int_{{\boldsymbol p},k}  & \frac{(2\pi)^2 \delta(k_0) \delta(p^2)}{((p_0^2 - ({\boldsymbol p} - {\boldsymbol k})^2 )^2 + \varepsilon^2) ( ({\boldsymbol k}{}^2)^2 + \varepsilon^2 )} \big( - 2 (d-1) p_0^3 \big) \big( p_0^2 - ({\boldsymbol p} - {\boldsymbol k})^2 \big)   \\
    & \times \bigg[ {\boldsymbol k}^2 + \frac12 \frac{\varepsilon \lambda^3 ({\boldsymbol k}{}^2)^4 }{\lambda^4 ({\boldsymbol k}{}^2)^4 + (1 - 2\lambda^2 {\boldsymbol k}{}^2 )^2 \varepsilon^2} -\frac12 \frac{\varepsilon^3 \lambda {\boldsymbol k}{}^2 (1-2\lambda^2 {\boldsymbol k}{}^2 ) }{\lambda^4 ({\boldsymbol k}{}^2)^4 + (1 - 2\lambda^2 {\boldsymbol k}{}^2 )^2 \varepsilon^2}  \bigg]\,, \nonumber
\end{align}
where we have assumed $p_0>0$ and taken $\frac{\varepsilon}{{(p^2)}^2 + \varepsilon^2} \to \pi \delta(p^2)$ from the start since this does not cause any singularity in the calculation.
The first term (i.e., the ${\boldsymbol k}^2$ term) in the square bracket gives the Feynman gauge result. The other two terms, together with the overall factor $( ({\boldsymbol k}{}^2)^2 + \varepsilon^2 )^{-1}$, give a contribution that is smaller than
$ {\boldsymbol k}{}^2 \delta({\boldsymbol k}^2)$ in the limit $\varepsilon\to0$, which is automatically vanishing. This can be seen by noting that
\begin{align}
    & \frac{1}{({\boldsymbol k}{}^2)^2 + \varepsilon^2} \frac{\varepsilon \lambda^2 ({\boldsymbol k}{}^2)^4 }{\lambda^4 ({\boldsymbol k}{}^2)^4 + (1 - 2\lambda^2 {\boldsymbol k}{}^2 )^2 \varepsilon^2} \leq \frac{\varepsilon \lambda^2 ({\boldsymbol k}{}^2)^2 }{\lambda^4 ({\boldsymbol k}{}^2)^4 + (1 - 2\lambda^2 {\boldsymbol k}{}^2 )^2 \varepsilon^2} \overset{\varepsilon \to 0}{\longrightarrow} \pi \lambda^2 ({\boldsymbol k}{}^2)^2 \delta( \lambda^2 ({\boldsymbol k}{}^2)^2 ) = \frac{\pi {\boldsymbol k}{}^2 }{2}  \delta({\boldsymbol k}{}^2) \nonumber \, , \\
    & \frac{1}{({\boldsymbol k}{}^2)^2 + \varepsilon^2} \frac{\varepsilon^3  {\boldsymbol k}{}^2 |1-2\lambda^2 {\boldsymbol k}{}^2 | }{\lambda^4 ({\boldsymbol k}{}^2)^4 + (1 - 2\lambda^2 {\boldsymbol k}{}^2 )^2 \varepsilon^2} \leq \frac{\varepsilon}{({\boldsymbol k}{}^2)^2 + \varepsilon^2} \frac{ {\boldsymbol k}{}^2  }{|1 - 2\lambda^2 {\boldsymbol k}{}^2|} \overset{\varepsilon \to 0}{\longrightarrow} \frac{\pi {\boldsymbol k}{}^2  }{|1 - 2\lambda^2 {\boldsymbol k}{}^2|} \delta({\boldsymbol k}{}^2) = \pi {\boldsymbol k}{}^2 \delta({\boldsymbol k}{}^2) \nonumber \, ,
\end{align}
and since the rest of the integral is non-singular at ${\boldsymbol k} = 0$ because the possible divergence in the principal value $\mathcal{P} \big[ (p_0^2 - ({\boldsymbol p} - {\boldsymbol k})^2 )^{-1} \big] $ is tamed by the ${\boldsymbol k}{}^2$ factor in the integral measure. Thus, at any finite $\lambda$, only the term that corresponds to the Feynman gauge calculation remains, which is still true in the limit $\lambda\to0$ if we take $\eta\to0$ first.

For completeness, we evaluate the remaining (finite) piece that gives the gauge invariant result for the difference. Now we can set $d=3$ because the integrals are strictly convergent and no regularization is required. The result is
\begin{align}
    g_E^{\rm Q\bar{Q}}(p_0) -  g_E^{\rm Q}(p_0) &= T_F g^4 N_c (N_c^2-1) \int_{{\boldsymbol p},k}  (2\pi)^2 \delta(k_0) \delta(p^2)  \mathcal{P} \left( \frac{ - 2 (d-1) p_0^3 }{(p_0^2 - ({\boldsymbol p} - {\boldsymbol k})^2 ) {\boldsymbol k}{}^2} \right) \nonumber \\
    &= \frac{T_F g^4 N_c (N_c^2-1)}{(2\pi)^3} \left( - 8 p_0^3 \right) \int_0^\infty \diff|{\boldsymbol p}| |{\boldsymbol p}|^2 \delta(p_0^2 - |{\boldsymbol p}|^2) \int_0^{\infty} \diff|{\boldsymbol k}| |{\boldsymbol k}|^2 \int_{-1}^1 \frac{\diff u}{ {\boldsymbol k}{}^2 (2 |{\boldsymbol p}| |{\boldsymbol k}| u - {\boldsymbol k}{}^2 ) } \nonumber \\
    &= \frac{T_F g^4 N_c (N_c^2-1) p_0^3}{(2\pi)^3} (-2) \int_0^\infty \frac{\diff |{\boldsymbol k}|}{|{\boldsymbol k}|} \ln \left| \frac{1 - |{\boldsymbol k}|/(2p_0)}{1+|{\boldsymbol k}|/(2p_0) } \right| \nonumber \\
    &= \frac{T_F g^4 N_c (N_c^2-1) p_0^3}{(2\pi)^3} 2 \int_0^\infty \frac{\diff x}{x} \ln \left| \frac{1 - x}{1+x } \right| \\\nonumber
    &= \frac{T_F g^4 N_c (N_c^2-1) p_0^3}{(2\pi)^3} \pi^2 \, ,
\end{align}
where we used our assumption $p_0>0$. Thus we reproduced the difference observed in Ref.~\cite{Binder:2021otw} by using mixed axial-Feynman gauge.

\end{widetext}

\end{document}