\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}