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


\section{Introduction \& Related Work}
\label{sec:introduction}

Logic synthesis is a crucial optimization step in the flow from the RTL description to the final GDS data. Minimizing the gate count of a logic circuit, even without timing constraints, is known to be NP-hard. With today's available compute power, logic functions with only a few inputs have come within reach of optimal synthesis (e.g., using SAT solvers). The synthesis of larger circuits, however, remains in the hand of algorithms relying on heuristics.

Deep learning has revolutionized the field of computer vision and many others by replacing engineered heuristics with multi-layered learned heuristics, leading to dramatic leaps in the quality of results (QoR). In the following, we investigate how the game of logic optimization can be taught to graph convolutional neural networks (GCNs) \cite{Kipf2017} operating on a graph representation of logic circuits using reinforcement learning (RL).

Two previous works have made first steps in this direction. DRiLLS \cite{Hosny2020} has applied RL at a coarse-grained level to select which commands/heuristics of UC Berkeley's ABC logic synthesis tool \cite{Epfl} should be applied, showing clear improvements in their QoR. Haaswijk \textit{et al.} \cite{Haaswijk2018} have investigated the direct application of RL to GCNs for logic optimization based on majority-inverter graphs (MIGs), showing that they can achieve optimality for small logic functions by comparing the results against synthesis based on SAT solvers. This method has the key limitations that
1) its embedding method does not scale to larger circuits, requiring a larger embedding dimension as the circuit size increases,
2) it has to be trained on a graph with an identical number of nodes as the one to be optimized, and
3) it only applies a single action at one node of the graph, resulting in $\mathcal{O}(N^2)$ compute time.

In this work, we combine RL with GCNs based on MIG logic graphs and focus on overcoming these limitations of \cite{Haaswijk2018}. Specifically, we contribute
1) a novel embedding method that is independent of the graph size while demonstrating its capability to achieve a high QoR;
2) a compact RL-based framework resulting in a GCN with only a few trainable parameters that enables inference on arbitrarily large graphs regardless of the graph size used during training;
3) a GCN that learns an optimized set of \textit{atomic} actions for each node, applying local modifications at many nodes within the graph with each inference step;
4) an evaluation of our method that includes 500-node random logic circuits.

\section{Proposed Method}
\label{sec:systemdesign}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/system_small.pdf}
\caption{Schema of the DRL system for logic optimization.}
\label{fig:system}
\end{figure}
Our method, shown in \cref{fig:system}, performs logic synthesis on circuits composed by 3-input majority and inverter gates represented in our model as MIGs, whose algebra allows defining a sound and complete move set \cite{Amaru2016}. At every time step, the system applies the instructions provided by an RL agent which learns autonomously which transformations are best to apply to each circuit configuration. The training process is organized by episodes with a fixed number of steps. At episode start, an initial graph is generated and a sequence with sets of transformations is applied until episode end. The RL agent learns its policy by interacting with the \textit{environment}. This is responsible for consistently applying the agent's chosen legal transformations while blocking illegal transformation attempts. The latter includes handling collisions of simultaneous actions from distinct nodes. To help the agent explore the space more efficiently, we provide additional inference rules which ease applicability, mitigating the reward sparsity problem. In this setup, it is no longer needed to commute each node inputs in specific ways to allow action application. The \textit{environment} computes and sends to the agent the reward (size difference between initial and final graphs) achieved in an episode. Thus, all actions taken at episodes with size reduction are positively reinforced.

For a given node $n_{i}$ of a MIG and a maximum adjacency depth $d_{adj}$, we define the neighborhood $N_{i}$ of $n_{i}$ as the set of nodes $n_{j}$ which are reachable from $n_{i}$ traversing no more than $d_{adj}$ edges. We define a state space $S$ where $s\in S$ represents any subgraph that can be observed by the GCN in $N_{i}$. We call \emph{action} any graph transformation that keeps intact the truth table of the MIG (legal move). These actions are selected by the agent ($\Omega$ actions), or are always applied after each step ($\Lambda$ actions), cf. \cref{tab:graphTransforms}. Additionally, the identity action $\Omega \cdot I$ is always applicable.

\begin{table}
    \centering
    \caption{Local Majority-Inverter Graph Transformations}
    \label{tab:graphTransforms}
    \resizebox{\linewidth}{!}{
    \begin{tabular}{@{}l@{\hspace{1mm}}l@{\hspace{1mm}}l@{}}
        \toprule
        $\Omega \cdot C$ & Commutativity & $M(x,y,z) = M(y,x,z) = M(z,y,x)$ \\
        $\Omega \cdot A$ & Associativity & $M(x,u,M(y,u,z) = M(z,u,M(y,u,x))$ \\
        $\Omega \cdot A_{c}$ & Compl. Assoc. & $M(x,u,M(y,u',z) = M(x,u,M(y,x,z))$ \\
        $\Omega \cdot D$ & Distributivity & $M(x,y,M(u,v,z) = M(M(x,y,z),M(x,y,v),z)$ \\
        $\Omega \cdot I$ & Inverter Prop. & $M'(x,y,z) = M(x',y',z')$ \\
        \midrule
        $\Lambda \cdot M$ & Majority & $M(x,x,z) = x$ and $M(x,x',z) = z$ \\
        $\Lambda \cdot R$ & Redundancy &
        $\forall u, v \in $ nodes: if in($u$)=in($v$), replace all $u$ with $v$\\
        \bottomrule
    \end{tabular}
    }
\end{table}

At every step of an episode, the agent observes the state of each node $n_{i}$ in the MIG by performing a GCN forward pass. At each hidden layer, the network collects for each $n_{i}$ features vectors (FVs) from adjacent nodes $n_{k}$. The FVs of its 3 inputs are concatenated with 3 aggregated FVs for the output. The nodes connected to the output are binned based on its input port (0, 1, or 2) of the connection, then the FVs in each of the bins are summed up to produce the aggregated FVs. The resulting FVs are passed through linear layers and non-linearities, while the process is repeated for every layer. In this layout, the network explores for each $n_{i}$ a broader $N_{i}$ at each layer. Finally, the network outputs a distribution with action probabilities for each node given its state $P(a\mid s)$, and the agent samples an action for each node. In accordance, the environment applies the legal selected moves and this process is repeated until the episode ends, when the reward is computed and used to update the model parameters. We train our model with Policy Gradient, as the sharp changes on estimation policy require the explicit state exploration from on-policy algorithms. They also break ties on equivalent sequences by exploiting one of them.

\begin{figure}
\centering
\includegraphics[width=0.97\linewidth]{./figures/nodeID_scalable_embeddings.pdf}
\caption{Representation of the pattern seen by the model within a node neighborhood for (1) traditional \textit{Node-ID} and (2) our proposed \textit{scalable} embeddings. For (1), the model explores neighborhoods by comparing feature vectors through paths in all directions, while the latter unravels sequences of edges through closed loops of length up to $d_{adj}$. Below, the unraveled gathered sequences of features are shown, which encode the inverters and gates traversed by the path. }
\label{fig:embeddings}
\end{figure}

\textbf{Model embeddings:} A major obstacle of GCN-based models is the need to embed each node with the necessary information. Common options are learned embeddings based on the node type, which lacks relevant information, and node IDs, which are non-scalable as they depend on the graph size.
We propose an alternative, scalable embedding method in which every node can just ``see'' whether a node in its neighborhood is ``itself'' or ``not itself''. This implies that every node neighborhood is initialized with its node-specific embedding before a forward pass. In the following, we demonstrate that this is sufficient to learn the required features for logic optimization. \cref{fig:embeddings} depicts the feature collection process.

The proposed method nevertheless has some limitations. 1) The multi-action schema makes it impractical to use actions involving an arbitrary number of nodes, such as substitution or relevance. These have been proven to be effective in heuristic optimization approaches. However, the system can still get to any state without these actions in a few more steps \cite{Amaru2016}. 2) The specific embeddings for each node have a memory and compute time footprint depending on the graph's connectivity. However, this is not a scalability issue as it does not depend on graph size.

\section{Results}

 We train our model with 3-input and 4-input decompositions obtained through sum-of-products (SOP) and compare the mean size reduction obtained with the maximum achievable, computed with \textit{cirkit}'s exact synthesis \cite{Soeken}. We also train with two sets of 1000 random graphs with 100 inputs, 2 outputs and sizes $(50, 500)$. Finally, we train our model with \textit{C1355} and \textit{C880} benchmarks \cite{Iscas85}. Size reductions achieved with Random graphs and benchmarks are compared with the ones obtained with ABC \cite{Epfl}. The following table shows the average size reduction achieved with each dataset and method, and their relative improvement with respect to the baseline.

\begin{table}[h!]
\label{tab:meansizered}
\centering

\resizebox{0.9\linewidth}{!}{
\begin{tabular}{ lrrrrr }
\toprule
Dataset & baseline$^a$ & \multicolumn{2}{c}{------ DLLOA \cite{Haaswijk2018} ------} & \multicolumn{2}{c}{------ \textbf{ours} ------} \\
 & MSR$^b$ & MSR$^b$ & relMSR$^c$ & MSR$^b$ & relMSR$^c$\\
\midrule
SOP$_{3}$  & 8.65 & \textbf{8.65} & \textbf{100\%} & 7.38 & 85\% \\
SOP$_{4}$  & 24.59 & 20.41 & 83\% & \textbf{24.57} & \textbf{100\%}\\
Rand 50 & 25.54 & -- & -- & \textbf{44.92} & \textbf{175\%} \\
Rand 500 & 271.24 & -- & -- & \textbf{413.68} & \textbf{152\%} \\
C880\cite{Iscas85}  & \textbf{39} & -- & -- & 17 & 44\%\\
C1355\cite{Iscas85}  & 114 & 98 & 86\% & \textbf{106} & \textbf{93\%}\\
\bottomrule
\end{tabular}
}

\vspace{1mm}
\raggedright
\footnotesize{$^a$ SOP datasets: optimal synthesis (cirkit \cite{Soeken}); all others: ABC \cite{Epfl}.\\
$^b$ Mean Size Reduction: Average difference between initial and final graphs.\\
$^c$ Relative Mean Size Reduction: Average size reduction relative to baseline.\\
}
\end{table}

The model provides competitive results for the benchmarks and SOP datasets while outperforming ABC for randomly generated graphs. This difference in performance is influenced by the fact that ABC operates with AIGs, which resembles the original structure of benchmarks and SOPs, while our LO operates with MIGs, used for generating random graphs.

\section{Conclusion}
\label{sec:conclusion}
We have proposed to combine GCNs with RL with a novel node embedding to perform logic optimization on MIGs. We show promising results on several benchmark graphs, outperforming ABC with a size reduction of up to 1.75$\times$ on larger, random graphs.

\bibliographystyle{IEEEtran}