



\subsection{Impact parameter space}

In order to solve for the splitting function in Eq.~(\ref{eq:AMY}), we switch from transverse momentum  $\mathbf{p}_{\perp}$ to transverse impact parameter $\bb_\perp$ space according to
\begin{eqnarray}
\mathbf{f}_{(z,P)}(\bb_\perp) = \frac{1}{16 P^2z^2(1{-}z)^2} \int \! \frac{\mathrm{d} ^2\mathbf{p}_{\perp}}{(2\pi)^2}~e^{i \bb_\perp \mathbf{p}_{\perp}}\mathbf{g}_{(z,P)}(\mathbf{p}_{\perp})
\end{eqnarray}
such that
\begin{eqnarray}
\mathbf{g}_{(z,P)}(\mathbf{p}_{\perp}) = 16 P^2z^2(1{-}z)^2~\int \! \mathrm{d} ^2\bb_\perp ~e^{-i \bb_\perp \mathbf{p}_{\perp}}~\mathbf{f}_{(z,P)}(\bb_\perp) 
\end{eqnarray}
where for later convenience we have absorbed a pre-factor $16 P^2 z^2 (1{-}z)^2$ into the definition, we need to calculate
\begin{eqnarray}
\frac{\mathrm{d} \Gamma_{ij}}{\mathrm{d} z} (P,z) =\alpha_s P_{ij}(z)\text{Im} \left[ 8\nabla_{\bb_\perp} \cdot \mathbf{f}_{(z,P)}(\bb_\perp) \right]
\end{eqnarray}
we get $\mathbf{p}_{\perp} \to -i \nabla_{\bb_\perp}$ and $\mathbf{p}_{\perp}^2 \to -  \nabla_{\bb_\perp}^2$ such that
\begin{eqnarray}
\frac{-2i}{16 P^2z^2(1{-}z)^2}~\nabla_{\bb_\perp} \delta^{(2)}(\bb_\perp) &=& i\left( M_{\rm eff}(z,P)- \frac{\nabla_{\bb_\perp}^2}{2 P z(1{-}z)} \right)~\mathbf{f}_{(z,P)}(\bb_\perp) \nonumber \\
&&+ \left[ C_{1} \bar{C}(\bb_\perp) + C_{z} \bar{C}(z\bb_\perp) + C_{1-z} \bar{C}((1{-}z)\bb_\perp) \right]~\mathbf{f}_{(z,P)}(\bb_\perp)\;, \nonumber \\\label{eq:ODE} 
\end{eqnarray}
where $\bar{C}(\mathbf{b})$ denotes the elastic scattering rate in impact parameter space stripped of its color factor.


\subsection{General strategy}
Since for $\bb_\perp \to 0$ the contributions from the collision term vanish, the structure of the solution near the origin is already contained in the "free" solution $\mathbf{f}^{(0)}_{(z,P)}(\bb_\perp)$. By matching the most singular terms near the origin one has
\begin{eqnarray}
\frac{-2i}{16 P^2z^2(1{-}z)^2} \nabla_{\bb_\perp} \delta^{(2)}(\bb_\perp) &=&- i\frac{\nabla_{\bb_\perp}^2}{2 P z(1{-}z)}~\mathbf{f}^{(0)}_{(z,P)}(\bb_\perp)\;,
\end{eqnarray}
which yields
\begin{eqnarray}
\lim_{\bb_\perp \to 0} \mathbf{f}^{(0)}_{(z,P)}(\bb_\perp) =  \frac{1}{8\pi P z(1{-}z)} \frac{\bb_\perp}{|\bb_\perp|^2} 
\end{eqnarray}
as can be easily seen from considering the two dimensional version of Gauss law. Denoting $b_\perp=|\bb_\perp|$ and expressing 
\begin{eqnarray}
\mathbf{f}^{(0)}_{(z,P)}(\bb_\perp) = \bb_\perp f^{(0)}_{(z,P)}(b_\perp)\;,
\end{eqnarray}
such that component wise
\begin{eqnarray}
\nabla_{\bb_\perp}^{i} \nabla_{\bb_\perp}^{i}  \bb_\perp^{j} f^{(0)}_{(z,P)}(b_\perp) = \bb_\perp^{j} \left( f^{(0)''}_{(z,P)}(b_\perp)+ \frac{3}{b_\perp} f^{(0)'}_{(z,P)}(b_\perp) \right)
\end{eqnarray}
we then have to solve
\begin{eqnarray}
\left[\partial_{b_\perp}^{2} + \frac{3}{b_\perp} \partial_{b_\perp} - 2P z (1{-}z) M_{\rm eff}(z,P) \right] f^{(0)}_{(z,P)}(b_\perp) = 0\;,
\end{eqnarray}
for $b_\perp>0$ with the boundary conditions
\begin{eqnarray}
\lim_{b_\perp \to 0}  {f}^{(0)}_{(z,P)}(b) = \frac{1}{8\pi P z (1{-}z) b_\perp^2} \;, \qquad  \lim_{b_\perp \to \infty}  {f}^{(0)}_{(z,P)}(b_\perp)=0\;,
\end{eqnarray}
 to determine the "free'' solution. By implementing the correct boundary conditions, one finds that the appropriate solution to Bessel's equation is given by
 \begin{equation}
 {f}^{(0)}_{(z,P)}(b_\perp)=  \frac{1}{8\pi P z (1{-}z) b_\perp^2} ~ \sqrt{2 P z(1{-}z) M(z,P)}b_\perp~K_{1}\Big(\sqrt{2 Pz(1{-}z)M(z,P)} b_\perp \Big)\;.
 \end{equation}
 where $K_{\nu}(z)$ denotes the modified Bessel function of rank $\nu$ and argument $z$. Note that the free solution $ {f}^{(0)}_{(z,P)}(b_\perp)$ is entirely real and hence does not contribute to the expression for the radiation rate.
 

Based on the "free'' solution ${f}^{(0)}_{(z,P)}(b_\perp)$, the solution $f_{(z,P)}(b_\perp)$ to the full evolution equation can be conveniently expressed as
 \begin{eqnarray}
 f_{(z,P)}(b_\perp)= {f}^{(0)}_{(z,P)}(b_\perp)+  {f}^{(1)}_{(z,P)}(b_\perp) \;,
 \end{eqnarray}
 where ${f}^{(1)}_{(z,P)}(b_\perp)$ describes the modifications due to the elastic scattering kernel $\bar{\Gamma}(b_\perp)$ and satisfies the linear inhomogeneous differential equation
 \begin{align}
 & \left[\partial_{b_\perp}^{2} + \frac{3}{b_\perp} \partial_{b_\perp} 
 - 2 P z(1{-}z) M_{\rm eff}(z,P) 
 \right. \nonumber \\ & \left. \phantom{\frac 12} {}
 +  2 i  P z(1{-}z) \Big[ C_{1} \bar{C}(b) + C_{z} \bar{C}(z b_\perp) + C_{1{-}z} \bar{C}((1{-}z)b_\perp) \Big]  \right]
 {f}^{(1)}_{(z,P)}(b) 
\nonumber \\ = & -2i  P z(1{-}z) \Big[ C_{1} \bar{C}(b_\perp) + C_{z} \bar{C}(z b_\perp) + C_{1-z} \bar{C}((1{-}z)b_\perp) \Big] {f}^{(0)}_{(z,P)}(b_\perp)\;.
 \end{align}
Since the singular part of $ f_{(z,P)}(b_\perp)$ is already captured by the contribution ${f}^{(0)}_{(z,P)}(b_\perp)$, one concludes that the function ${f}^{(1)}_{(z,P)}(b_\perp)$ has to be regular in the vicinity of the origin, i.e.
\begin{eqnarray}
\lim_{b_\perp \to 0} {f}^{(1)}_{(z,P)}(b_\perp) = c_{1} \,.
\end{eqnarray}
Now as we have already seen from the free case, the above differential equation clearly features diverging solutions in the limit $b_\perp \to 0$. Hence the requirement that ${f}^{(1)}_{(z,P)}(b_\perp)$ remains finite fixes one of the integration constants in the general solution of the differential equation. Similarly, the second integration constant can be determined from the requirement that the solution remains regular in the limit $b_\perp \to \infty$, i.e.
\begin{eqnarray}
\lim_{b_\perp \to \infty} {f}^{(1)}_{(z,P)}(b_\perp) =0\;.
\end{eqnarray}
Expressing the general solution of the linear inhomogeneous ODE as the general solution to the homogeneous ODE plus a special solution to the inhomogenous ODE and choosing the special solution to be regular for $b \to 0$, we can formally write
\begin{eqnarray}
{f}^{(1)}_{(z,P)}(b_\perp)= c_{\rm div} {f}^{(1) {\rm hom, div}}_{(z,P)}(b_\perp) + c_{\rm reg} {f}^{(1) {\rm hom, reg}}_{(z,P)}(b_\perp) +{f}^{(1) {\rm inhom, reg}}_{(z,P)}(b_\perp)
\end{eqnarray}
where the first coefficient $c_{\rm div}$ multiplying the divergent contribution  ${f}^{(1) {\rm hom, div}}_{(z,P)}(b_\perp)$ is fixed as $c_{\rm div}=0$ from the requirement that ${f}^{(1)}_{(z,P)}(b_\perp)$ is regular. Conversely the second coefficient $c_{\rm reg}$ can be fixed from the requirement that ${f}^{(1)}_{(z,P)}(b_\perp)$ remains finite in the limit $b_\perp \to \infty$, i.e.
\begin{eqnarray}
\label{eq:cReg}
 c_{\rm reg} = - \lim_{b_\perp \to \infty} \frac{{f}^{(1) {\rm inhom, reg}}_{(z,P)}(b_\perp)}{{f}^{(1) {\rm hom, reg}}_{(z,P)}(b_\perp) }
\end{eqnarray}
Hence, if we choose to solve the differential equations with the following initial conditions in the limit $b\to0$
\begin{align}
{f}^{(1) {\rm hom, reg}}_{(z,P)}(b_\perp=0) & =g_{\mathrm{s}}^2 T\;, & {f'}^{(1) {\rm hom, reg}}_{(z,P)}(b_\perp=0) & = 0 \;,
\\
{f}^{(1) {\rm inhom, reg}}_{(z,P)}(b=0) & = g_{\mathrm{s}}^2 T\;, &
{f'}^{(1) {\rm inhom, reg}}_{(z,P)}(b_\perp=0) & = 0 \;,
\end{align}
we can exploit the fact that ${f}^{(1)}_{(z,P)}(b)$ and its derivatives do not diverge in the limit $b \to 0$ to determine the relevant derivative at the origin ($\nabla_{\bb_\perp} \bb_\perp = 2$) as
\begin{eqnarray}
 \text{Im}\left[8\nabla_{\bb_\perp} \cdot \mathbf{f}_{(z,P)}(\bb_\perp) \right]_{\bb_\perp=0} = 16 \text{Im}~c_{1} = 16~g_{\mathrm{s}}^2 T~\text{Im}~c_{\rm reg}\;,
\end{eqnarray}
with $ c_{\rm reg}$ given by Eq.~(\ref{eq:cReg}). By combining all the relevant pieces the splitting rate is then determined by
\begin{eqnarray}
\frac{\mathrm{d} \Gamma_{ij}}{\mathrm{d} z} (P,z) =\frac{4}{\pi}~g_{\mathrm{s}}^4 T P_{ij}(z)~\text{Im}~c_{\rm reg} \,.
\end{eqnarray}

\subsection{Numerical implementation}
We solve the homogeneous and inhomogeneous ODE for ${f}^{(1) {\rm hom, reg}}$ and ${f}^{(1) {\rm inhom, reg}}$ numerically using an adaptive fourth order Runge-Kutta scheme. We start at a finite value of $b_\perp = b_\perp^{\rm min}$, where the initial conditions for ${f}^{(1) {\rm hom, reg}}$ are determined from the series expansion of the solution around $b_\perp=0$ as
\begin{eqnarray}
{f}^{(1) {\rm hom, reg}}(b_\perp^{\rm min}) &=& g_{\mathrm{s}}^2 T \left( 1 + \frac{1}{4} M_{\rm eff} P z \bar{z} b_\perp^{\rm min\; 2} \right) \;,  
\end{eqnarray}
and similarly for the inhomogeneous solution ${f}^{(1) {\rm inhom, reg}}$. Since the latter is sensitive to the elastic scattering kernel, different expansions have to be employed for the different interaction kernels. Specifically, for the LO potential, the small $b_\perp$ expansion takes the form
\begin{align}
    {f}^{(1) {\rm inhom, reg}}(b_\perp^{\rm min}) & - {f}^{(1) {\rm hom, reg}}(b_\perp^{\rm min}) = 
    \frac{-i c(T)}{32\pi}(g_{\mathrm{s}} T b_\perp^{\rm min})^2\ln(g_{\mathrm{s}} T b_\perp^{\rm min}) \left( C_1 {+} z^2 C_z {+} \bar{z}^2 C_{\bar{z}} \right) \nonumber\\
    &-\frac{i(g_{\mathrm{s}} T b_\perp^{\rm min})^2}{128\pi} \left( 
        \left[4 \frac{\hat{q}_0}{4 g_{\mathrm{s}}^4 T^3} + c(T) (4\ln(2)- 3) \right] (C_1+z^2 C_z + \bar{z}^2 C_{\bar{z}}) \right.
        \nonumber\\
        & \left. \phantom{ \frac{i(g_{\mathrm{s}} T b_\perp^{\rm min})^2}{128\pi}} 
        \hspace{1.5em} {}+ 4c(T)\left(C_1 \ln\tfrac{1}{2} +z^2 C_z \ln\tfrac{z}{2} + \bar{z}^2 C_{\bar{z}}\ln\tfrac{\bar{z}}{2}\right)
     \right)\;,
\end{align}
where $ c(T) = - \frac{C_\mathrm{R}}{8 \pi} \frac{\zeta(3)}{\zeta(2)} \left( -\frac{1}{2 g_{\mathrm{s}}^2} + \frac{3 y}{2} \right) (g_{3\mathrm{d}}^2 b_\perp)^2 $.
Based on the terms in the expansion, the starting value $b_\perp^{\rm min}$ is then chosen as
\begin{equation}
b_\perp^{\rm min}=10^{-4} \text{min}\left(\frac{1}{\sqrt{ M_{\rm eff} P z \bar{z} }},\frac{1}{8 \pi  P z \bar{z}} ,  (g_{\mathrm{s}} T)^{-1}, \frac{1}{\sqrt{(C_{1} +C_{z} z^2+ C_{\bar{z}} \bar{z}^2) (g_{\mathrm{s}} T)^{2}}} \right)\;,
\end{equation}
such that higher order corrections are effectively negligible at double precision level. Starting from the initial conditions at $b=b_{\rm min}$, we then solve the differential equations for ${f}^{(1) {\rm hom, reg}}$ and ${f}^{(1) {\rm inhom, reg}}$, up to a maximal value of $b_\perp=b_\perp^{\rm max}$, which is chosen similarly according to the criterion that
\begin{equation}
|f^{(0)}(b_\perp^{\rm max})/{f}^{(1) {\rm inhom, reg}}(b_\perp^{\rm max})| < 10^{-16}\,.
\end{equation}
Based on the solution, we then determine the constant $c_{\rm reg}$ from the ratio 
\begin{equation}
 c_{\rm reg} = - \frac{{f}^{(1) {\rm inhom, reg}}_{(z,P)}(b_\perp^{\rm max})}{{f}^{(1) {\rm hom, reg}}_{(z,P)}(b_\perp^{\rm max})}\;.
\end{equation}
Such that the differential radiation rate is ultimately given by
\begin{equation}
\frac{\mathrm{d} \Gamma_{ij}}{\mathrm{d} z} (P,z) =\frac{4}{\pi}~g_{\mathrm{s}}^4 T P_{ij}(z)~\text{Im}~c_{\rm reg}\;.
\end{equation}



\subsection{Effective splitting rate }
The starting point for the calculation of the inelastic splitting rates  $\mathrm{d} \Gamma_{ij}/ \mathrm{d} z(P,z)$ is the following expression
\begin{eqnarray}
\frac{\mathrm{d} \Gamma_{ij}}{\mathrm{d} z} (P,z) = \frac{\alpha_s P_{ij}(z)}{[2Pz(1{-}z)]^2} \int \! \frac{\mathrm{d} ^2 \mathbf{p}_{\perp}}{(2\pi)^2} ~\text{Re} \left[ 2\mathbf{p}_{\perp} \cdot \mathbf{g}_{(z,P)}(\mathbf{p}_{\perp}) \right]\;,
\end{eqnarray}
where we follow the notation of P.~Arnold in Appendix A of \cite{Arnold:2008iy}. Here $P_{ab}(z)$ are the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) splitting functions
\begin{eqnarray}
    P_{gg}(z) &=& 2 C_\mathrm{A} \frac{[1-z(1{-}z)]^2}{z(1{-}z)}\;,\quad
    P_{qg}(z)  = C_\mathrm{F} \frac{1+(1{-}z)^2}{z}\;,\quad
    P_{gq}(z) = \frac{1}{2} \left(z^2+(1{-}z)^2\right) \;. \nonumber\\
\end{eqnarray}
The function $\mathbf{g}_{(z,P)}(\mathbf{p}_{\perp})$, which encodes the current-current correlator, satisfies the following integral equation:
\begin{eqnarray}
    \label{eq:AMY}
    2\mathbf{p}_{\perp} &=& i \delta E(z,P,\mathbf{p}_{\perp}) \mathbf{g}_{(z,P)}(\mathbf{p}_{\perp}) + \int \! \frac{\mathrm{d} ^2\mathbf{q}_\perp}{(2\pi)^2}~\bar{C}(q_\perp)   \\
    && \hspace{120pt} \times \left\{ C_{1} \left[ \mathbf{g}_{(z,P)}(\mathbf{p}_{\perp}) - \mathbf{g}_{(z,P)}(\mathbf{p}_{\perp} -\mathbf{q}_\perp) \right] \right. \nonumber \\
    && \hspace{140pt} \left. + \, C_{z} \left[ \mathbf{g}_{(z,P)}(\mathbf{p}_{\perp}) - \mathbf{g}_{(z,P)}(\mathbf{p}_{\perp} -z\mathbf{q}_\perp) \right] \right. \nonumber \\
    && \hspace{140pt} \left. + \, C_{1-z} \left[ \mathbf{g}_{(z,P)}(\mathbf{p}_{\perp}) - \mathbf{g}_{(z,P)}(\mathbf{p}_{\perp} -(1{-}z)\mathbf{q}_\perp) \right] \right\} \, . \nonumber
\end{eqnarray}
The energy difference $\delta E(z,P,\mathbf{p}_{\perp})$ is written
\begin{eqnarray}
\delta E(z,P,\mathbf{p}_{\perp}) = \frac{\mathbf{p}_{\perp}^2}{2Pz(1{-}z)} +M_{\rm eff}(z,P)\;,
\end{eqnarray}
where $M_{\rm eff}(z,P)$ is given in terms of the asymptotic masses $m^2_{\infty,(1,z,1-z)}$ of the particles with momentum fractions $1,z,1-z$ as
\begin{eqnarray}
M_{\rm eff}(z,P)=  \frac{m^2_{\infty,(z)}}{2zP}  + \frac{m^2_{\infty,(1{-}z)}}{2(1{-}z)P} -\frac{m^2_{\infty,(1)}}{2P}.
\end{eqnarray}
For the asymptotic masses we use the leading order results given by
\begin{eqnarray}
    m^2_{\infty,g}=\frac{m_{\mathrm{D}}^2}{2}= \frac{g_{\mathrm{s}}^2T^2}{6}\left( C_\mathrm{A} +\frac{N_{\mathrm{f}}}{2}\right)\;,\qquad
    m^2_{\infty,q}= C_\mathrm{F} \frac{g_{\mathrm{s}}^2 T^2}{4}\;.
\end{eqnarray} 
The color factors are given by
\begin{eqnarray}
\label{colorfactors}
C_{1}=\frac{1}{2} \left( C^{R}_{z} + C^{R}_{1-z} - C^{R}_{1} \right)\;, \nonumber &\quad&
C_{z}=\frac{1}{2} \left( C^{R}_{1-z} + C^{R}_{1} - C^{R}_{z} \right)\;,  \nonumber \\
C_{1-z}&=&\frac{1}{2} \left( C^{R}_{1} + C^{R}_{z} - C^{R}_{1-z} \right)\;, 
\end{eqnarray}
where $C^{R}_{(1,z,1-z)}$ denote the Casimir of the representation of the particle carrying momentum fraction $1,z,1-z$, i.e. $C^{R}=C_\mathrm{F}$ for quarks and $C^{R}=C_\mathrm{A}$ for gluons. Since the color factors have been factored out, the rate $\bar{C}(q)$ in Eq.~(\ref{eq:AMY}) denotes the elastic scattering rate stripped of its color factor.%
\footnote{%
Technically, we make an approximation when we use $\bar{C}(q_\perp)$
and \eqref{colorfactors} in \eqref{eq:AMY}.  Our expression
for $\bar{C}(q_\perp)$ is for a 
fundamental-representation particle picking up
transverse momentum $q_\perp$ in the medium, and is determined
by the Wilson line we described;
but to compute the evolution of a mixed state
with an emitter in one amplitude and an emitter and a gluon
in the conjugate amplitude, we really need a 3-Wilson-line
object, as described in Ref.~\cite{CaronHuot:2008ni}.
This factorizes into a combination of Wilson-line-pair
contributions as shown in \eqref{colorfactors} up to corrections
which are at most NNLO.  Theoretically, this three-Wilson-line
object could also be directly computed in EQCD, but it
would be more challenging, and in particular we would
need a separate numerical calculation for each pair $(b_\perp,x)$.
Here we are assuming that the NNLO corrections are small
and the 3-Wilson-line object factorizes as shown.
As far as we know, \textsl{every} treatment in the literature
makes this same approximation!
}

We solve Eq.~(\ref{eq:AMY}) in impact-parameter space for the splitting rates using the numerical procedure outlined in the Appendix~\ref{sec:numerical} to obtain the rates of medium induced $g \to gg$, $q \to qg$ and $g \to q\bar{q}$ splittings and supply the software as part of the arXiv submission package. Since the fully resummed AMY rate includes both the high-energy limit where the Landau-Pomeranchuk-Migdal (LPM) effect is prominent, and the low-energy limit which follows a Bethe-Heitler (BH) rate, we will briefly discuss these limits before addressing the numerical results in more detail.


\subsection{Bethe-Heitler regime}
 When the typical momentum of splitting is small $(Pz(1{-}z) \ll \omega_{\rm BH} \sim T )$, the formation time of the radiation is small and interference between scatterings can be neglected. In this so-called Bethe-Heitler regime, one can then solve the rate Eq.~(\ref{eq:AMY}) in an opacity expansion, corresponding to expansion in the number of elastic scatterings with the medium. By consider the limit of a single scattering, we obtain the following semi-analytic expressions for the rates (c.f. Appendix \ref{sec:BH})
 \begin{align}
    \label{eq:BHRate}
    \frac{\mathrm{d} \Gamma_{g\to gg}^{BH}}{\mathrm{d} z} (P,z) = & 
    g_{\mathrm{s}}^4  T P_{gg}(z) \times \\ 
    & \left[  \frac{C_A}{2} Q\left( \tilde{\mu}^{2}_{g\to gg}(z) \right)  + \frac{C_\mathrm{A}}{2} Q\left( \tilde{\mu}^{2}_{g\to gg}(z)/z^2 \right) + \frac{C_\mathrm{A}}{2} Q\left( \tilde{\mu}^{2}_{g\to gg}(z)/\bar{z}^2 \right) \right], \nonumber \\ 
    \frac{\mathrm{d} \Gamma_{q\to gq}^{BH}}{dz} (P,z) = & 
    g_{\mathrm{s}}^4 T P_{qg}(z) \times \nonumber \\
    & \left[\frac{C_\mathrm{A}}{2} Q\left( \tilde{\mu}^{2}_{q\to gq}(z) \right)  +(C_\mathrm{F}-\frac{C_\mathrm{A}}{2}) Q\left( \tilde{\mu}^{2}_{q\to gq}(z)/z^2 \right) +\frac{C_\mathrm{A}}{2} Q\left( \tilde{\mu}^{2}_{q\to gq}(z)/\bar{z}^2 \right) \right], \nonumber \\ 
    \frac{\mathrm{d} \Gamma_{g\to qq}^{BH}}{\mathrm{d} z} (P,z) = &
     g_{\mathrm{s}}^4 T  P_{gq}(z) \times \nonumber \\
     & \left[ (C_\mathrm{F} - \frac{C_\mathrm{A}}{2}) Q\left( \tilde{\mu}^{2}_{g\to qq}(z) \right) + \frac{C_\mathrm{A}}{2} Q\left( \tilde{\mu}^{2}_{g\to gq}(z)/z^2 \right) +\frac{C_\mathrm{A}}{2} Q\left( \tilde{\mu}^{2}_{q\to gq}(z)/\bar{z}^2 \right) \right],\nonumber
\end{align}
where, denoting $a=m^2_{\infty,q}/m_D^2$, one has
\begin{align}
\tilde{\mu}^{2}_{g\to gg}(z)&=\frac{1-z(1{-}z)}{2}\;,  & \tilde{\mu}^{2}_{q\to gq}(z)&=\frac{z}{2}+a(1{-}z)^2\;, & \tilde{\mu}^{2}_{g\to qq}(z)& =\frac{2a-z(1{-}z)}{2}\;,
\end{align}
and
\begin{eqnarray}
\label{eq:QBHDef}
Q(\tilde{\mu}^2)=  \frac{m_{\mathrm{D}}^2}{2\pi g_{\mathrm{s}}^2 T} \int \! \frac{\mathrm{d} ^2 \mathbf{p}_{\perp}}{(2\pi)^2} \int \! \frac{\mathrm{d} ^2\mathbf{q}_\perp}{(2\pi)^2} ~ \bar{C}(m_D\mathbf{q}_\perp)  \left[ \frac{\mathbf{p}_{\perp}}{\mathbf{p}_{\perp}^2+\tilde{\mu}^2} - \frac{(\mathbf{p}_{\perp}-\mathbf{q}_\perp) }{(\mathbf{p}_{\perp}-\mathbf{q}_\perp)^2+\tilde{\mu}^2}  \right]^2\;.
\end{eqnarray}
While the above relation is formulated in momentum space, the  integral defining $Q(\tilde{\mu}^2)$ in Eq.~(\ref{eq:QBHDef}) can also be evaluated using the kernel in position space as show in Appendix \ref{sec:BH}.


\subsection{Deep LPM regime}
Conversely, in the limit of a very high-energy parton $(P\gg T)$ traversing a thick medium, the typical number of rescatterings within the formation time of bremsstrahlung can be large, indicating that interferences between many soft scatterings which contribute to the total transverse momentum transfer during the formation of the radiation need to be considered. Simplifications occur in the limit $Pz(1{-}z) \gg \omega_{\rm BH} \sim T$, where the splitting probes the small $b_\perp$ behavior of the momentum broadening kernel which can be expressed as
\begin{eqnarray}
C(b_\perp) = -\frac{g_{\mathrm{s}}^4 T^3}{16\pi}\mathcal{N}b_\perp^2 \log( \xi m_{\mathrm{D}}^2 b_\perp^2/4)
\end{eqnarray}
where $\mathcal{N} = \frac{\zeta (3)}{\zeta(2)} \left( 1 + \tfrac{N_{\mathrm{f}}}{4} \right)$. In accordance with the discussion in Sec.~\ref{sec:nonpert_Cbp}, the coefficient $\frac{g_{\mathrm{s}}^4T^2}{16\pi}\mathcal{N} b_\perp^2$ gives the leading logarithmic behavior $b_\perp \log(b_\perp^2)$ and the coefficient $\xi$ captures the $b_\perp^2$ behavior. Specifically at for the LO kernel $\xi_{\rm LO}=e^{2\gamma_\mathrm{E}-2}\simeq 0.429313$ can be determined analytically, while for the NLO and non-perturbative kernels, we obtain $\xi_{\rm NLO} \simeq 1.355 \cdot10^{-3}$ and $\xi_{\rm NP} = 4\frac{g_{\mathrm{s}}^4 T^2 }{m_{\mathrm{D}}^2} e^{-4\pi\frac{  \hat{q}_0}{g_{\mathrm{s}}^4 T^3 \mathcal{N}}}\simeq 0.1702$ from a fit of the small $b_\perp$ behavior. Following \cite{Arnold:2008zu}, the rate equation can be solved iteratively in an inverse logarithmic expansion to obtain\footnotemark
\footnotetext{ The rate in our notation is related to the notation in \cite{Arnold:2008zu} as $\frac{\mathrm{d} \Gamma_{a\to bc}}{\mathrm{d} z}(P) = \frac{(2\pi)^3}{P\nu_a}\gamma^{a}_{bc}(P|zP,\bar{z}P)$ with $\nu_g = 2(N_\mathrm{c}^2-1)$ and $\nu_q = 2N_\mathrm{c}$.}
\begin{eqnarray}
    \label{eq:LPMRate}
    \frac{\mathrm{d} \Gamma_{a\to bc}}{\mathrm{d} z}(P)&=&\frac{g_{\mathrm{s}}^2}{16\pi^2\sqrt{2} Pz(1-z)}~P_{ab}(z)~m_{\mathrm{D}}^2 \mu^2_\perp(P,z)\;,
\end{eqnarray}
where $\mu^2_\perp(P,z)$ is self-consistently determined from 
\begin{eqnarray}
    \mu^2_\perp(P,z) &=& \frac{g_{\mathrm{s}}^2 T^2\mathcal{N}}{m_{\mathrm{D}}^2}\frac{g_{\mathrm{s}} T}{m_{\mathrm{D}}} \left( \frac{2}{\pi} z(1{-}z)\frac{P}{T} \right)^{1/2}
    \bigg( C_1 \log\left( \frac{\alpha  \mu_\perp^2}{\xi}\right)  \nonumber \\
    &&\qquad \left. + \, C_z z^2 \log\left(\frac{\alpha  \mu_\perp^2}{\xi z^2}\right) + \, C_{1-z} (1-z)^2 \log\left(\frac{\alpha  \mu_\perp^2}{\xi(1{-}z)^2}\right)\right)^{1/2}\;,
\end{eqnarray}
with  $\alpha= e^{\gamma_e+\pi/4}$.


\subsection{Results}


\begin{figure}[h!]
    \includegraphics[width=0.98\textwidth]{Plots/RateCompare-T500.pdf}
    \caption{Splitting rate at $T=250$MeV (dashed blue lines) and $T=500$MeV (full purple lines) for the processes $g \rightarrow gg$ (top),  $q \rightarrow gq$ (middle), $g \rightarrow qq$ (bottom). Different columns correspond to parent energies $p=10T$ (left) and $p=100T$ (right). We compare with rates computed using the perturbative leading order (orange) and next-to-leading order (green) elastic broadening kernels. The Bethe Heitler rates and LPM rates are shown with dashed lines and circles respectively, using the color of the corresponding kernel. }
    \label{fig:SplittingRates500}
\end{figure}

\begin{figure}[h!]
    \includegraphics[width=0.98\textwidth]{Plots/RateMomentum.pdf}
    \caption{Momentum dependence of the splitting rate at $T=500$MeV (dashed blue lines) for the processes $g \rightarrow gg$ (top),  $q \rightarrow gq$ (right), $g \rightarrow qq$ (bottom). Dashed lines and open circles correspond to the approximate rates in Bethe Heitler regime (\ref{eq:BHRate}) and the deep LPM regime (\ref{eq:LPMRate}). }
    \label{fig:MomentumRate}
\end{figure}

Numerical results for the medium-induced splitting rates are presented in Fig.~\ref{fig:SplittingRates500}, where we show the rates for the non-perturbative broadening kernel $C_{\rm QCD}(b_\perp)$ at $T=250,500$MeV (left/right columns), along with the corresponding results obtained for the leading order ($C_{\rm QCD}^{\rm LO}(b_\perp)$) and next-to-leading order ($C_{\rm QCD}^{\rm NLO}(b_\perp)$) determinations of the collisional broadening kernel.
Different panels in Fig.~\ref{fig:SplittingRates500} show the rates for the $g\to gg$ (top),$q \to qg$ (center) and $g \to q\bar{q}$ (bottom) processes, for different parton energies $p=10T$ (left) and $p=100T$ (right).  The momentum dependence of the rate is shown in Fig.~\ref{fig:MomentumRate}, for the non-perturbative kernel at $T=500$MeV, using the color palette in a logarithmic scale to distinguish between different momentum of the parent particle $p=10-1000T$. In both figures, we also show the Bethe-Heitler rates in Eq.~(\ref{eq:BHRate}) (dashed lines) and the deep LPM rates (circles) Eq.~(\ref{eq:LPMRate}).

Starting with the rates in Fig.~\ref{fig:SplittingRates500}, one observes that the non-perturbatively determined splitting rates for the two different temperatures do not display any remarkable difference, leading to basically the same emission rates in units of $[g^4T]$. As expected, the momentum dependence of the rate clearly displays LPM suppression at large typical momentum $Pz(1{-}z) \gg T$ as well as an unsuppressed Bethe-Heitler rate in the other limit $Pz(1{-}z) \ll T$ seen in Fig.~\ref{fig:MomentumRate}.
We also observe that at low energy $z(1{-}z)E \ll T$ where the large impact parameter (small momentum transfer) is more important, the non-perturbative result is closer to the NLO rate as they both have a similar behavior at large impact parameter. 
Despite the apparent non-convergence of the perturbative series for $C(\bbp)$ observed in \cite{CaronHuot:2008ni}, it is interesting to observe that the splitting rates for the non-pertubatively determined $C(\bbp)$ mostly fall between the LO and NLO results and it would be interesting to compare them to higher order perturbative calculations in the future.
Conversely, in the LPM suppressed regime at high energy  $z(1{-}z)E \gg T$, where small impact parameter (large momentum transfer) is relevant, the non-perturbative rate is closer to the LO rate which again agrees with the behavior of the elastic kernel (cf. 
Fig.~\ref{fig:Spline}). 


\section{Introduction}
\label{sec:intro}
\input{Introduction.tex}

\section{Non-perturbative broadening kernel}
\label{sec:nonpert_Cbp}
\input{nonpert_Cbp.tex}

\section{Medium induced splitting rates} 
\label{sec:inelastic_rates}
\input{Splitting.tex}


\section{Conclusion}
\label{sec:conclusion}
\input{Conclusion.tex}

\section*{Acknowledgements}
We would like to express our gratitude to Jacopo Ghiglieri for fruitful discussions in the early stages of this problem,
for careful comments on the original draft,
and for providing us with a Mathematica notebook performing the numerical Fourier transform of \eqref{Cqp_hard}.
We thank Shuzhe Shi for comments on the draft.
This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the CRC-TR 211 ’Strong-interaction matter under extreme conditions’– project number 315477589 – TRR 211. N.\ S.\ acknowledges support from Academy of Finland grants 267286 and 320123.


\subsection{Matching for the transverse momentum broadening kernel}
Since EQCD is a low-energy effective theory for full QCD, we know that the two of them should agree%
\footnote{We use the terms IR and UV to refer to small and large momentum scales; for coordinate-space regimes, we refer to short-distance and long-distance.}
in the infrared (IR) regime, well below the hard scale $p \ll \pi T$. However, as one goes further into the ultra-violet (UV), discrepancies should arise.  Therefore, the EQCD result $C_\mathrm{EQCD}(\bbp)$ determined in~\cite{Moore:2020wvy} cannot agree with its full QCD counterpart $C_\mathrm{QCD}(\bbp)$, which is our ultimate object of desire.
In order to ensure the correct short-distance behavior and to keep the fully nonperturbative long-distance information from lattice EQCD at the same time, we schematically express
\begin{equation}	\label{match_strategy1}
C_\mathrm{QCD}(\bbp) = \left( C_\mathrm{QCD}(\bbp) - C_\mathrm{EQCD}(\bbp) \right) + C_\mathrm{EQCD}(\bbp) \, .
\end{equation}
The quantity in parenthesis is then computed as part
of a matching calculation between full QCD and EQCD.
We will carry out this matching calculation in transverse
momentum space and Fourier transform the result to $b_\perp$-space.
Such matching calculations are free of large IR effects
because the IR regime by definition agrees between the
two theories.  Therefore the relevant scale in the
matching calculation is $2\pi T$, where the effective
coupling is (barely) perturbative and perturbation
theory should be applicable. Also, the effects of large
statistical functions, which make IR perturbative behavior
poorly behaved in thermal QCD even at weak coupling,
do not arise in this matching calculation.
Therefore we can determine this quantity perturbatively.
For the last EQCD term, we use the nonperturbative
lattice results:
\begin{equation}
\label{eq:match2}
C_{\mathrm{QCD}}(b_\perp) \approx \left( C_{\mathrm{QCD}}^{\mathrm{pert}}(b_\perp) - C_{\mathrm{EQCD}}^{\mathrm{pert}}(b_\perp) \right) + C_{\mathrm{EQCD}}^{\mathrm{latt}}(b_\perp) \,.
\end{equation}


Calculations of the perturbative terms in Eq.~(\ref{eq:match2}) are currently available at $\mathcal{O}(g_{\mathrm{s}}^4)$,%
\footnote{Note that $C_{\mathrm{EQCD}}^{\mathrm{pert}}(b_\perp)$ and $C_{\mathrm{QCD}}^{\mathrm{pert}}(b_\perp)$ obey slightly different power counting schemes; while the expansion parameter in the former case is the (four-dimensional) strong coupling constant $g_{\mathrm{s}}^2$, the latter case features the (three-dimensional) expansion parameter $g_{3\mathrm{d}}^2$. Up to subleading corrections, these two expansion parameters are related via $g_{3\mathrm{d}}^2 \approx g_{\mathrm{s}}^2 T$. The full-QCD scale-dependence enters EQCD through the perturbative matching procedure \cite{Kajantie:1997tt,Laine:2005ai} that relates the EQCD effective parameters to the full-QCD coupling $g_{\mathrm{s}}$ at a renormalization scale $\Lambda_{\overline{\mathrm{MS}}} \approx 341~\mathrm{MeV}$ \cite{Bruno:2017gxd}.} 
but only in transverse momentum $q_\perp$ space.
Therefore we will write the expression in parenthesis
in \eqref{eq:match2} in $\mathbf{q}_\perp$ space and then perform
a Fourier transform to $\bb_\perp$ space.
The perturbative QCD contribution for $q_\perp \gg m_D$ is given by \cite{Arnold:2008vd}
\begin{equation}	\label{Cqp_hard}
C^{\mathcal{O}(g_{\mathrm{s}}^4)}_\mathrm{QCD}(q_\perp) = \frac{g_{\mathrm{s}}^4 T^3 C_\mathrm{R}}{q_\perp^4} \!\int \! \frac{\mathrm{d}^3 p}{(2 \pi)^3} \frac{p - p_z}{p} 
\left[ 2 C_\mathrm{A} n_\mathrm{B}(p) \left( 1 {+} n_\mathrm{B}(p') \right) + 4 \, N_{\mathrm{f}} T_{\mathrm{f}} \, n_\mathrm{F}(p) \left( 1 {-} n_\mathrm{F}(p') \right) \right] .
\end{equation}
Here $C_\mathrm{R}$ is the Casimir of the jet constituent's representation $R$ of $SU(3)$, $C_\mathrm{A}=3$ is the respective Casimir operator of the adjoint representation, $T_\mathrm{f}=1/2$ is the normalization of the fundamental representation of $SU(3)$, $N_{\mathrm{f}}$ is the number of massless fermion flavors,
and $p$ and $p' = p + \frac{\boldsymbol q_\perp^2 + 2 \boldsymbol q_\perp \hspace{-1pt} \cdot \, {\boldsymbol p}}{2 (p - p_z)}$ are the momenta of the medium particle before/after undergoing the scattering. While \eqref{Cqp_hard} provides the correct UV limit of $C(q_\perp)$ in QCD, the (unphysical) IR limit of \eqref{Cqp_hard} coincides with the (unphysical) UV limit of EQCD, both from perturbation theory \cite{CaronHuot:2008ni} and the lattice \cite{Moore:2019lgw}.
We quote the full perturbative $C_\mathrm{EQCD}(\bbp)$ at leading order \cite{Aurenche:2002wq} and next-to-leading order \cite{CaronHuot:2008ni} in \eqref{CEQCDbp_NLO} of Appendix~\ref{sec:pert_rate}. 
Here, however, we only need the large-$q_\perp$
expansion to NLO, that is, to order $q_\perp^{-3}$,
since this order matches the precision of the full-QCD
calculation \cite{Ghiglieri:2018ltw}.  To this order,
the EQCD value reads
\begin{equation}	\label{Cqp_subtr}
C^\mathrm{pert}_\mathrm{EQCD}(q_\perp) \xrightarrow{q_\perp \gg m_{\mathrm{D}}} C^\mathrm{pert}_\mathrm{subtr}(q_\perp) = \frac{C_\mathrm{R} g_{\mathrm{s}}^2 T m_{\mathrm{D}}^2}{q_\perp^4} - \frac{C_\mathrm{R} C_\mathrm{A} g_{\mathrm{s}}^4 T^2}{16 \, q_\perp^3} \, .
\end{equation}
The first term in \eqref{Cqp_subtr} cancels against the IR limit of \eqref{Cqp_hard}, avoiding double-counting degrees of freedom in marrying the soft and hard contributions to $C_\mathrm{QCD}(\bbp)$.  The second term turns into a positive
linear term in $C(\bbp)$.  This removes the negative
linear small-$\bb_\perp$ behavior found in EQCD, ensuring
the positivity of the full $C(\bbp)$.
This feature was first pointed out by
Caron-Huot \cite{CaronHuot:2008ni}.

The lattice data $C_{\mathrm{EQCD}}^{\mathrm{latt}}$
of~\cite{Moore:2020wvy} is obtained in position $(b_\perp)$ space and a direct numerical Fourier transform to momentum space is a delicate issue.  Therefore we choose to
perform the matching calculation in $\bb_\perp$ space,
meaning that we need to Fourier transform the $q_\perp$-space
expressions we have just presented
using \eqref{subtraction_FT}, and perform the matching entirely in position space. Schematically, rewriting \eqref{match_strategy1} into
\begin{equation}	\label{match_strategy2}
C_{\mathrm{QCD}}(b_\perp) \approx \int \! \frac{\mathrm{d}^2 q_\perp}{(2\pi)^2} \left( C_{\mathrm{QCD}}^{\mathrm{pert}}(q_\perp) - C^{\mathrm{pert}}_{\mathrm{subtr}}(q_\perp) \right)  \left( 1 - e^{i\boldsymbol q_\perp \cdot \mathbf{b}_{\bot}} \right) + C_{\mathrm{EQCD}}^{\mathrm{latt}}(b_\perp) 
\end{equation}
is more accurate, where the Fourier transform of the difference of \eqref{Cqp_hard} and \eqref{Cqp_subtr} requires a numerical integration but is feasible with standard tools, for instance Mathematica 12. 

Even though only the two lowest temperatures, $T = 250\, \mathrm{MeV}$ and $T = 500\,\mathrm{MeV}$, will be directly relevant for the subsequent computation of radiative emission rates, we provide the fully matched results for $C_{\mathrm{QCD}}(b_\perp)$ in Tab.~\ref{tab:full_Cbp} for all temperatures for which $C_{\mathrm{EQCD}}^{\mathrm{latt}}(b_\perp)$ was reported in~\cite{Moore:2019lgw}. Our results are cast into dimensionless ratios with the help of the three-dimensional coupling $g_{3\mathrm{d}}^2$, and plotted in Fig.~\ref{fig:asymptotics}, where we present $C_\mathrm{QCD}(\bbp)$ at the two temperatures of further relevance for this work together with the limiting infrared and ultraviolet behavior that we will discuss further below.

\begin{table}[ptbh] 
\centering {\small
\begin{tabular}{|c|c|c|c|c|}	
\hline
 & & & & \\[-10pt]
$g_{3\mathrm{d}}^2 b_\perp$ & $\left. \frac{C(\bbp) }{ g_{3\mathrm{d}}^2} \right\vert^{N_\mathrm f = 3}_{250~\mathrm{MeV}}$ & $\left. \frac{C(\bbp) }{ g_{3\mathrm{d}}^2} \right\vert^{N_\mathrm f = 3}_{500~\mathrm{MeV}}$ & $\left. \frac{C(\bbp) }{ g_{3\mathrm{d}}^2} \right\vert^{N_\mathrm f = 4}_{1~\mathrm{GeV}}$ & $\left. \frac{C(\bbp) }{ g_{3\mathrm{d}}^2} \right\vert^{N_\mathrm f = 5}_{100~\mathrm{GeV}}$ \\
 & & & & \\[-10pt]
\hline
$0.125$ & $-0.0011(44)$ & - & - & $-0.001(19)$ \\
$0.25$ & $0.0000(36)$ & $-0.0041(36)$ & $-0.0006(34)$ & $0.011(32)$ \\
$0.5$ & $0.00552(63)$ & $0.01244(87)$ & $0.02198(58)$ & $0.04553(49)$ \\
$0.75$ & $0.0242(11)$ & $0.0357(17)$ & $0.0559(10)$ & $0.09742(81)$ \\
$1.0$ & $0.03685(82)$ & $0.06181(61)$ & $0.09118(36)$ & $0.14850(33)$ \\
$1.5$ & $0.1025(17)$ & $0.1440(11)$ & $0.19449(94)$ & $0.28356(62)$ \\
$2.0$ & $0.1783(33)$ & $0.2467(28)$ & $0.3148(18)$ & $0.4304(12)$ \\
$2.5$ & $0.2747(41)$ & $0.3542(51)$ & $0.4443(32)$ & $0.5868(22)$ \\
$3.0$ & $0.3784(86)$ & $0.4756(46)$ & $0.5781(67)$ & $0.7425(45)$ \\
$4.0$ & $0.514(42)$ & $0.711(33)$ & $0.853(28)$ & $1.095(17)$ \\
$5.0$ & $0.815(94)$ & $1.017(96)$ & $1.09(10)$ & $1.359(13)$ \\
$6.0$ & $1.31(11)$ & $1.44(10)$ & $1.67(15)$ & $1.748(36)$ \\
\hline
$g^2$ & $3.725027$ & $2.763516$ & $2.210169$ & $1.066561$ \\
$y$ & $0.452423$ & $0.586204$ & $0.823449$ & $1.64668$ \\
\hline
$\hat{q}_{0} / g_{3\mathrm{d}}^6$ & $0.1465(78)$ & $0.185(10)$ & $0.3136(60)$ & $0.5665(47)$ \\
\hline
\end{tabular} }
\caption{
  Results for $\frac{C_\mathrm{QCD}}{g_{3\mathrm{d}}^2}(b_\perp)$ for four 
  temperatures and a range of transverse separations $b_\perp$.  
}
\label{tab:full_Cbp}
\end{table}

When comparing the lattice EQCD result $C_\mathrm{EQCD}^\mathrm{latt}(b_\perp)$ in Tab.~2 of \cite{Moore:2019lgw} and the fully matched QCD result $C_\mathrm{QCD}(\bbp)$ in our Tab.~\ref{tab:full_Cbp}, one finds that the matching only introduces marginal corrections at large $b_\perp$.
However, our matching procedure in \eqref{match_strategy2} largely cures the negative dip in the short-distance regime of EQCD, where its impact is much more significant than at large distances, as expected.
All but one small-$b_\perp$ values are consistent with $0$ within a single standard deviation, one value is a bit more than one standard deviation away from the positive region.
Nevertheless, not all central values of $C(\bbp)$ are manifestly positive.  This poses a numerical challenge later on,
and it also reflects the growing difficulty of extracting
precise data from EQCD at these short distances.
At such short distances we are instead better off trying
to handle EQCD perturbatively, so we will switch at
small $b_\perp$ to the asymptotic small-separation
perturbative expression in EQCD.
Similarly, the lattice data become noisy at very
large separation, and it is best to fit them to the
expected analytical large-$b_\perp$ form and to use this
analytical expression at the largest separations.
In these two limiting cases, the functional forms are available and are discussed in the following two subsections. We note that, whenever we speak of a long-distance limit in the following, we mean the limit of $b_\perp \gg \xi_\mathrm{max}$, where $\xi_\mathrm{max}$ is the longest correlation length in the theory, typically $\xi_\mathrm{max} \sim \frac{1}{g_{3\mathrm{d}}^2}$.
Similarly, whenever we speak of the short-distance limit, we mean $b_\perp \ll \xi_\mathrm{min}$, typically $\xi_\mathrm{min} \sim \frac{1}{m_{\mathrm{D}}}$ the smallest correlation length, such that the short-distance limit also genuinely surpasses the hard scale $b_\perp < (\pi T)^{-1}$, meaning that we exceed the range of validity of EQCD and probe the matching part in this limit.


\begin{figure}[htbp!] 
\centering
	 \subfigure[Short-distance limit]{\includegraphics[scale=0.55]%
	 {Plots/small_r_asymptotics.pdf}} 
	 \subfigure[Long-distance limit]{\includegraphics[scale=0.55]%
	 {Plots/large_r_asymptotics.pdf}} 
	\caption{Short-distance (a) and long-distance (b) behavior of  the momentum broadening kernel $C_\mathrm{QCD}(b_\perp)$. 
		Dashed lines correspond the asymptotic functional forms in \eqref{IR_limit} and \eqref{UV_limit}}
	\label{fig:asymptotics}
\end{figure}



\subsection{Long-distance limit of $C(\bbp)$}
We first focus on the long-distance limit of $C_\mathrm{QCD}(\bbp)$, which turns out to be the less complex issue. Beyond the known area-law form of $C_\mathrm{EQCD}(\bbp)$, subleading asymptotic corrections are found to be important to match to the numerical data at the values of $b_\perp$, where lattice calculations are feasible. Therefore, the infrared limit of the full $C(\bbp)$ reads
\cite{Laine:2012ht}
\begin{equation}	\label{IR_limit}
\frac{C_\mathrm{QCD}}{g_{3\mathrm{d}}^2} (b_\perp) \xrightarrow{b_\perp \gg \; 1/g_{3\mathrm{d}}^2} A + \frac{\sigma_\mathrm{EQCD}}{g_{3\mathrm{d}}^4} g_{3\mathrm{d}}^2 b_\perp +  \frac{g_{\mathrm{s}}^4 C_\mathrm{R}}{\pi} \left[ \frac{y}{4} \left( \frac{1}{6} - \frac{1}{\pi^2} \right) + \frac{C_\mathrm{A}}{8 \pi^2 g_{\mathrm{s}}^2} \right] \log (g_{3\mathrm{d}}^2 b_\perp) \, ,
\end{equation}
where $\sigma_\mathrm{EQCD}$ is the string tension of EQCD \cite{Laine:2005ai} and $A$ is a scale-setting fitting constant. The (subleading) $\log (\dots)$-part originates from the modified Fourier-transform of the $\frac{1}{q_\perp^2}$-term in the small-$q_\perp$ expansion of \eqref{Cqp_hard}, already provided in \cite{Arnold:2008vd}. Since it is introduced by the matching procedure, this term does not show up in the IR limit of $C_\mathrm{EQCD}(\bbp)$ \cite{Moore:2019lgw}. Meanwhile, we have re-expressed all occurrences of $N_{\mathrm{f}}$ in \eqref{IR_limit} in terms of the 1-loop expression of $m_{\mathrm{D}}^2 = \frac{(2 C_\mathrm{A} + N_{\mathrm{f}})}{6} g^2 T^2$ and $y = \frac{m_{\mathrm{D}}^2}{g_{3\mathrm{d}}^4} \Big|_{\bar{\mu}=g_{3\mathrm{d}}^2}$. For $y$, in turn, we used the two-loop expression from \cite{Laine:2005ai}, giving numerical values in Tab.~\ref{tab:full_Cbp}. This procedure does not spoil the perturbative power-counting since the difference is formally subleading. In practice, we find better agreement with our data using this procedure which can be thought of as a selective resummation of higher-order contributions to the screening mass $m_{\mathrm{D}}^2$. 

We plot the limiting behavior in \eqref{IR_limit} against our data for the two lowest temperatures $T = 250,\,500\,\mathrm{MeV}$ in Fig.~\ref{fig:asymptotics}(b). As one can see from Fig.~\ref{fig:asymptotics}(b), the linear contribution to \eqref{IR_limit} dominates numerically already in the displayed window of $g_{3\mathrm{d}}^2 b_\perp$. However, we find the logarithmic contribution, though numerically small, has still a non-negligible impact on our result. Just as in the case of EQCD, we see the onset of the asymptotic behavior at smaller $g_{3\mathrm{d}}^2 b_\perp$ the larger the temperature becomes, i.e.\ the smaller the (running) four-dimensional coupling $g_{\mathrm{s}}$ is.


\subsection{Short-distance limit of $C(\bbp)$}
Evaluating $C_\mathrm{EQCD}^\mathrm{latt}(b_\perp)$ at smaller $b_\perp$ demands smaller lattice spacings $a$, as discretization errors of the operator occur in powers of $\frac{a}{b_\perp}$~\cite{DOnofrio:2014mld}. As a consequence, the signal-to-noise ratio shrinks due to critical slowing down as one approaches $b_\perp \hspace{-3pt} \to 0$. Additionally, the leading term in the small-$b_\perp$-expansion of $C_\mathrm{EQCD}(\bbp)$ is precisely an EQCD artifact and therefore cancelled by the matching, diminishing the physically relevant information further without reducing the noise. Thus, the urge for an effective analytical description is clear.

The short-distance limit of $C_\mathrm{QCD}(\bbp)$ precisely gives rise to the momentum diffusion coefficient $\hat{q}$
\begin{equation}
\frac{C_\mathrm{QCD}}{g_{3\mathrm{d}}^2} (b_\perp) \xrightarrow{b_\perp \ll \; 1/m_{\mathrm{D}}} 
\frac{1}{4} \frac{\hat{q}}{g_{3\mathrm{d}}^6}(g_{3\mathrm{d}}^2 b_\perp)^2 \, ,
\end{equation}
where $\hat{q}$ consists of a scale-dependent logarithm and constant part. Clearly, at very small $b_\perp$, the logarithm is supposed to be the dominant part, and can be extracted by a modified Fourier-transform of the UV limit of \eqref{Cqp_hard} as described below. However, in keeping only the leading logarithmic terms, we do not find a smooth connection to our data, even at the smallest-$b_\perp$ data points. There are two ways to resolve this issue: Either computing data points at even smaller $b_\perp$ from the lattice or improving on the analytical side of the limiting behavior. Option one is certainly possible but numerically extremely costly due to the requirement of finer lattices combined with critical slowing down. Option two is somewhat easier to realize and will therefore be further pursued. 

By taking into account the constant part of the momentum diffusion coefficient, $\hat{q}_0$, which receives both perturbative and non-perturbative contributions, the short-distance limit of $C_\mathrm{QCD}(\bbp)$ reads
\begin{equation}	\label{UV_limit}
\frac{C_\mathrm{QCD}}{g_{3\mathrm{d}}^2} (b_\perp) \xrightarrow{b_\perp \ll \; 1/m_{\mathrm{D}}} - \frac{C_\mathrm{R}}{8 \pi} \frac{\zeta(3)}{\zeta(2)} \left( -\frac{1}{2 g_{\mathrm{s}}^2} + \frac{3 y}{2} \right) (g_{3\mathrm{d}}^2 b_\perp)^2 \log(g_{3\mathrm{d}}^2 b_\perp) + \frac{1}{4} \frac{\hat{q}_0}{g_{3\mathrm{d}}^6} (g_{3\mathrm{d}}^2 b_\perp)^2 \, ,
\end{equation}
where we applied the resummation of subleading contributions to the screening mass $m_{\mathrm{D}}^2$ just as in the IR case.

$\hat{q}_0$ would be most straight-forwardly determined from a fit to a few of the smallest-$b_\perp$-data-points of $C_\mathrm{QCD}(\bbp)$. However, the cancellations between the EQCD result and the matching contribution make a naive fit numerically quite unstable and badly-constrained. Also, the lattice EQCD results at small $b_\perp$ are strongly correlated among each other, which translates to $C_\mathrm{QCD}(\bbp)$. 
Determining $\hat{q}_0$ directly from the lattice-EQCD-data and its correction from $\left(C_\mathrm{QCD}^{\rm pert}-C_\mathrm{subtr}^{\rm pert}\right)(b_\perp)$ separately turns out to be a numerically more robust procedure, leading to the results quoted in Tab.~\ref{tab:full_Cbp}. We note that the momentum broadening coefficient of EQCD has already been determined from a small-$b_\perp$-fit in \cite{Moore:2019lgw}, also taking into account the correlation between the data points. We find that, together with the matching contribution, the values of $\hat{q}_0$ in Tab.~\ref{tab:full_Cbp} are smaller than the lattice-EQCD results of~\cite{Moore:2019lgw}.

We present a comparison of our result for $C_\mathrm{QCD}(b_\perp)$ in Tab.~\ref{tab:full_Cbp} to the short-distance asymptotics \eqref{UV_limit} in Fig.~\ref{fig:asymptotics}(a). Since in the accessible range of $b_\perp$, deviations between the asymptotic behavior and the data points are still sizeable, the question where to switch from our data points to the asymptotic description is in fact non-trivial and introduces an uncertainty into the calculation of radiative emission rates. When parametrizing $C_\mathrm{QCD}(b_\perp)$, we will therefore  keep the switching point as a parameter and finally investigate its impact on the broadening kernel.


\subsection{Numerical interpolation of lattice data}
\input{Spline.tex}


