\section{Introduction}

The production of energetic hadrons and jets in heavy ion collisions is markedly different from the production of energetic electroweak bosons. The latter clearly exhibit ``binary scaling'': weak bosons and high-energy photons are produced as if nucleons from each nucleus were independently undergoing inelastic binary collisions.
Hadron and jet measurements, on the other hand, display evident deviations from binary scaling.
These deviations are understood to be
a consequence of the formation of a quark-gluon plasma in relativistic nuclear collisions: energetic parton \emph{production} does follow ``binary scaling''; it is their subsequent interactions with the plasma that leads to parton energy loss, and consequently to an apparent deviation from binary scaling for hadronic observables.

This characteristic phenomena of jet and hadron ``energy loss'' in heavy ion collisions has been observed at both the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC)~\cite{Majumder:2010qh, Muller:2012zq, Mehtar-Tani:2013pia,Roland:2014jsa, Connors:2017ptx, Busza:2018rrf}. 
Energetic partons are produced at the earliest stage of heavy ion collisions, and they propagate through all the different phases of the collisions. 
As a consequence of their interactions
with the quark-gluon plasma, the momentum distribution of these energetic partons changes distinctly compared to the baseline observed in proton-proton collisions.\footnote{The role of hadronic energy loss is still under investigation. See Ref.~\cite{Dorau:2019ozd} for example.}
This makes them important probes of the properties of the deconfined nuclear plasma produced in heavy ion collisions.

A number of different formalisms have been used to model the interaction of energetic partons with the constituents of the plasma~\cite{Majumder:2010qh}. 
Fundamentally, most parton energy loss formalisms have a well-understood common core, yet applications to heavy ion collisions tend to require approximations and practical considerations that lead to non-negligible differences between parton energy loss models~\cite{Bass:2008rv,CaronHuot:2010bp,Armesto:2011ht}.
One difference between the models is the treatment of the underlying plasma, which
is often assumed
to be made of a large number of quarks and gluons with energies of $\lesssim 1$~GeV  in near local thermal equilibrium.
Whether these quarks and gluons are treated as dynamical entities or as static scattering centers is one of many differences in the energy loss formalisms~\cite{Bass:2008rv,CaronHuot:2010bp,Armesto:2011ht}.
The above assumption is important, given that  the quark-gluon plasma produced in heavy ion collisions is understood to be strongly coupled, and a  quasiparticle description may not be justified.\footnote{In particular,  hydrodynamic simulations of this plasma's evolution do not rely on a quasiparticle picture of deconfined nuclear matter until hadronization.}

A different phrasing of the above challenge is that even the energy loss of very energetic partons can be affected by non-perturbative effects from the strongly-coupled plasma.
Hard interactions --- those with large momentum transfer between the energetic parton and the plasma --- are expected to have smaller non-perturbative effects, or even be accessible perturbatively, as a consequence of the running of the QCD coupling.
On the other hand, ``soft'' parton-plasma interactions with small momentum transfer are expected to suffer the largest non-perturbative effects.\footnote{Note other works such as Refs.~\cite{Casalderrey-Solana:2014bpa,Casalderrey-Solana:2016jvj} assume that neither soft or hard interactions are perturbative, and consequently evaluate parton energy loss using gauge-field duality.}


A stochastic treatment of these soft interactions of energetic partons provides an alternative approach to account for non-perturbative effects --- an approach that is agnostic to the strongly- or weakly-coupled nature of the underlying deconfined plasma.
The dynamical details of the large number of soft interactions are encoded in a small number of transport coefficients.
The latter can be parametrized and constrained from measurements. They can also be studied using lattice techniques~(see for example Ref.~\cite[Section 4]{Ghiglieri:2015zma} and Refs.~\cite{Moore:2019lgw,Schlusser:2020lme} and references therein).
From a practical point of view, a stochastic description of a large number of soft interactions can also be more efficient numerically than a rate-based approach. 

A systematic hard-soft factorization of parton energy loss was proposed recently to describe parton propagation in a weakly-coupled QGP~\cite{Ghiglieri:2015ala,Ghiglieri:2015zma}. In this factorization, soft interactions are described as a stochastic process with drag and diffusion transport coefficients calculated perturbatively; hard interactions are solved with rates that are also calculated perturbatively. 
In the weakly-coupled regime, parton energy loss in this hard-soft factorizated scheme was shown to be equivalent to a fully rate-based treatment of parton energy loss~\cite{Ghiglieri:2015ala}. 
Importantly parton energy loss in this hard-soft factorization can also be extended to next-to-leading order~\cite{Ghiglieri:2015ala}, a feature beyond the scope of this work which we shall explore in the future.

The above work is based on the ``effective kinetic theory'' approach~\cite{Arnold:2002zm} derived for a weakly-coupled quark-gluon plasma.
In a weakly-coupled plasma, quark and gluon excitations are described as quasi-particles with effective properties related to the local density of the plasma. In this effective kinetic approach, the dynamics of quasi-particles are described by Boltzmann transport equations.
Leading order [$\mathcal{O}(\alpha_s)$] realizations of this effective kinetic approach --- extrapolated to large values of strong coupling constant $\alpha_s$ --- have been used widely to study parton energy loss~(see e.g. Refs.~\cite{Jeon:2003gi,Qin:2007rn, Schenke:2009gb,Burke:2013yra, Cao:2017zih}). 
Importantly, parton energy loss derived in this weakly-coupled effective kinetic theory can be related to other energy loss models that do not assume a weakly-coupled quark-gluon plasma.
Thus, while the current work is based on this weakly-coupled effective kinetic theory approach, there exists solid literature (see for example Ref.~\cite{CaronHuot:2010bp} and references therein) that can help relax  weak-coupling assumptions in the future.

In this work, we present the first numerical implementation of the hard-soft factorized parton energy loss model~\cite{Ghiglieri:2015ala} discussed above. 
For our implementation we utilize the public JETSCAPE framework~\cite{Putschke:2019yrg}, as it allows us a straightforward integration of our parton energy loss model with the other ingredients necessary for a full simulation of jet production in heavy ion collisions. 
We first test and validate this factorization of parton energy loss in the weak coupling regime for a static medium.

We introduce a dimensionless scale to quantify the kinematic range for which soft interactions can be described accurately with a stochastic approach.
We use this scale to discuss a hard-soft factorization model for a strongly-coupled quark-gluon plasma, relevant for phenomenological applications in heavy ion collisions.

Finally, we present an application of our new factorized model of parton energy loss by calculating the energy and fermion-number cascade of a hard parton propagating in a static medium, finding good agreement with known analytical approximations. 

\section{Hard-soft factorization of parton energy loss in the weakly-coupled regime: theory}

\subsection{Effective kinetic approach in weakly-coupled regime}

The evolution of an energetic parton in a thermal medium of temperature T can be described by a Boltzmann transport equation:
\begin{equation}
    \label{eq:Boltzmann}
    \left( \frac{\partial}{\partial t}+ \frac{\textbf{p}}{p} \cdot \mathcal{\nabla} \right) \delta f_a = - \mathcal{C}[\delta f_a, n_a]
   
   
\end{equation}
where $P=(p,\textbf{p})$ is the four-momentum of the energetic parton 
and $\mathcal{C}$ is the collision kernel of the parton with the medium.
We use the same notation for the parton momentum distributions as in Ref.~\cite{Ghiglieri:2015ala}:  the distribution of rare hard partons is $\delta f^a(\mathbf{p}, \mathbf{x}, t)$, to distinguish it from the quasi-thermal distribution of soft particles $n^a(\mathbf{p}, T(\mathbf{x}, t), \mathbf{u}(\mathbf{x}, t))$, where $\mathbf{u}$ is the flow velocity. In this notation, the total phase space distribution of quasiparticle $a$ is $f^a(\mathbf{p}, \mathbf{x}, T) = n^a(\mathbf{p}, T(\mathbf{x}, t), \mathbf{u}(\mathbf{x}, t))+\delta f^a(\mathbf{p}, \mathbf{x}, t)$.
We assume $p \gg T$ and $g \equiv \sqrt{\frac{\alpha_s}{4\pi}} \ll 1$. 
Because interactions between energetic partons themselves are rare and can be neglected, the Boltzmann equation is effectively linear in $\delta f_a$. 

At leading order, the interactions between quasi-particles can be divided as $2\leftrightarrow 2$ elastic interactions and $1\leftrightarrow 2$ inelastic interactions. Elastic $2\leftrightarrow 2$ processes refer to elementary scatterings involving two incoming particles and two outgoing particles without any radiation. Multiple soft $2\leftrightarrow 2$ scatterings between the energetic parton and the plasma can induce a collinear radiation. In the effective kinetic approach, these multiple soft scatterings are resummed consistently, to account for interference between subsequent collisions which lead to the Landau-Pomeranchuk-Migdal (LPM) effect.
This resummed collinear radiation is known as the effective $1\leftrightarrow 2$ process. 
The collision kernel of both $2\leftrightarrow 2$ and $1\leftrightarrow 2$ processes can be written as 
\begin{equation}
    \mathcal{C} = \mathcal{C}^{1\leftrightarrow 2}+\mathcal{C}^{2\leftrightarrow 2}. 
\end{equation}

Importantly, in our approach, we only follow the evolution of energetic partons with an energy above a cutoff $p_{\textrm{cut}} = 2$~GeV. 
Our assumption is that we should focus our efforts on high-$p_T$ observables which are dominated by partons above this cutoff.
After neglecting terms suppressed by $\exp(-p/T)$, the collision kernels $\mathcal{C}^{1\leftrightarrow 2}$ and $\mathcal{C}^{2\leftrightarrow 2}$ read: 
\begin{widetext}
\begin{equation}
\label{eq:1to2}
\begin{split}
    \mathcal{C}_a^{1\leftrightarrow 2}\left[\delta f\right] &= \frac{(2\pi)^3}{2|\bm{p}|^2\nu_a}\sum_{bc}\int^\infty_0dp'dq'\gamma_{bc}^{a}(\bm{p}; p'\hat{\bm{p}}, q'\hat{\bm{p}}')\delta(|\bm{p}|-p'-q')\\
    &\times \left\{\delta f^a(\bm{p})\left[1\pm n^b(p')\pm n^c(q')\right]-\left[\delta f^b(p'\hat{\bm{p}})n^c(q')+n^b(p')\delta f^c(k'\hat{\bm{p}})\right]\right\}\\
    &+\frac{(2\pi)^3}{|\bm{p}|^2v_a}\sum_{bc}\int^\infty_0dqdp'\gamma_{ab}^{c}(p'
    \hat{\bm{p}}; \bm{p}, q\hat{\bm{p}})\delta(|\bm{p}|+q-p')\\
    &\times \left\{\delta f^a(\bm{p})n^b(q)-\delta f^c(p'\hat{\bm{p}})\left[1\pm n^b(q)\right]\right\}\, , 
\end{split}
\end{equation}

\begin{equation}
\label{eq:2to2}
\begin{split}
    \mathcal{C}_a^{2\leftrightarrow 2}\left[\delta f\right] =& \frac{1}{4|\bm{p}|\nu_a}\sum_{bcd}\int_{\bm{k}\bm{p}'\bm{k}'}|\mathcal{M}_{cd}^{ab}(\bm{p}, \bm{k}; \bm{p}', \bm{k}')|^2(2\pi)^4\delta^{(4)}(P+K-P'-K')\\
    &\times \left\{\delta f^a(\bm{p})n^b(k)\left[1\pm n^c(p')\pm n^d(k')\right]-\delta f^c(\bm{p}')n^d(k')\left[1\pm n^b(k)\right]-n^c(p')\delta f^d(\bm{k}')\left[1\pm n^b(k)\right]\right\}\, , 
\end{split}
\end{equation}
\end{widetext}
where the notation for the Lorentz-invariant integration is
\begin{equation}
    \int_{\bm{k}} \dots \equiv \int \frac{d^3k}{2k(2\pi)^3}\dots
\end{equation}
and $\nu_a$ is the degeneracy of particle $a$. 

For $\mathcal{C}^{1\leftrightarrow 2}$, $a$ is the incoming hard parton with the momentum $\bm{p}$, and $b, c$ are outgoing particles with the momentum $\bm{p}', \bm{k}'$. $\gamma_{bc}^{a}$ is the splitting kernel of $a\rightarrow bc$, which can be calculated with the AMY integral equations~\cite{Arnold:2002zm,Jeon:2003gi}. 

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{plots/phase_division_elas.pdf}
    \includegraphics[width=\columnwidth]{plots/phase_division_inel.pdf}
    \caption{Treatment of different processes in the hard-soft factorized parton energy loss model}
    \label{fig:reformulation} 


\end{figure}

For $\mathcal{C}^{2\leftrightarrow 2}$, particle $a$ is the incoming energetic parton with momentum $\bm{p}$, particle $b$ is the plasma particle with the momentum $\bm{k}$ interacting with $a$, and particle $c, d$ are the outgoing particles with momentum $\bm{p}', \bm{k}'$. $\mathcal{M}_{cd}^{ab}$ is the matrix element of the elementary process $ab\rightarrow cd$. 

\subsection{Reformulating parton energy loss with hard-soft factorization}
\label{sec:factorized_model}

In the hard-soft factorized parton energy loss model introduced in Ref.~\cite{Ghiglieri:2015ala}, $1\leftrightarrow 2$ and $2\leftrightarrow 2$ processes are further divided into soft interactions and hard interactions. The collision kernel is rewritten as 
\begin{equation}
\mathcal{C} = \mathcal{C}^{1\leftrightarrow 2}_{\textrm{hard}}+\mathcal{C}^{1\leftrightarrow 2}_{\textrm{soft}}+\mathcal{C}^{2\leftrightarrow 2}_{\textrm{hard}}+\mathcal{C}^{2\leftrightarrow 2}_{\textrm{soft}}\, . 
\end{equation}

In this hard-soft factorized model, soft interactions described by $\mathcal{C}^{1\leftrightarrow 2}_{\textrm{soft}}$ and $\mathcal{C}^{2\leftrightarrow 2}_{\textrm{soft}}$ are treated stochastically with the Langevin equation. 

Hard inelastic interactions,  $\mathcal{C}^{1\leftrightarrow 2}_{\textrm{hard}}$, are treated with an emission rate as calculated from the AMY integral equations~\cite{Arnold:2002zm}.\footnote{We thank Guy D. Moore for his numerical solver for AMY integral equations.} We refer to them
as large-$\omega$ interactions.

The hard $2\leftrightarrow 2$ part $\mathcal{C}^{2\leftrightarrow 2}_{\textrm{hard}}$ is further divided as (i) large-angle interactions, and (ii) splitting approximation processes, based on the energy transfer $\omega$:  
\begin{equation}
\mathcal{C}^{2\leftrightarrow 2}_{\textrm{hard}} = \mathcal{C}_{\textrm{large-angle}}^{2\leftrightarrow2}+\mathcal{C}_{\textrm{split}}^{2\leftrightarrow2}\, . 
\end{equation}
The physical meaning of this separation is the following. Elastic collisions occur between an energetic parton ($p\gg T$) and a lower energy quasi-thermal quark or gluon ($k \sim T$). On rare occasions, the momentum transfer in these elastic collisions is sufficient to make the low-energy quark or gluon become an energetic parton with $k\gg T$; such partons are referred to in the literature as ``recoil partons''. The process through which a recoil parton is produced  is akin to a splitting process: a single energetic particle splits in two energetic ones. The kinematic of this process simplifies and it benefits from being treated separately.

The factorization of the phase space for this reformulation is summarized in Fig. \ref{fig:reformulation}. An ensemble of cutoffs are used to divide the different regions of phase space. We discuss the details of the different treatments and these cutoffs in the following subsections. 

\subsubsection{Treatment of hard interactions: inelastic case ($1\leftrightarrow 2$)}

\label{sec:factorized_model:inelastic_hard}

\begin{figure}
    \centering
    \includegraphics[width=0.8\linewidth]{plots/inel_diagram.pdf}
    \caption{Multiple soft scatterings induce the radiation of a soft gluon with energy $\omega$. We denote the radiations with $\omega > \mu_\omega$ as large-$\omega$ inelastic interactions.
   
    }
    \label{fig:inel_diagram}
\end{figure}

The diagram of a typical $1\leftrightarrow 2$ inelastic interaction is shown as Fig. \ref{fig:inel_diagram}. We assume the energy of the radiated particle is $\omega$. We define a hard-soft cutoff $\mu_\omega$ based on the radiated energy $\omega$, to divide $\mathcal{C}^{1\leftrightarrow 2}_{\textrm{hard}}$ and $\mathcal{C}^{1\leftrightarrow 2}_{\textrm{soft}}$. 
In the weakly-coupled regime, the cutoff $\mu_\omega$ is limited to $\mu_\omega \lesssim T$, where $T$ is the temperature of the thermal medium. 

Collinear radiations with energy $\omega  > \mu_\omega $ are included into the hard part, $\mathcal{C}^{1\leftrightarrow 2}_{\textrm{hard}}$; 
they are treated as usual with emission rates calculated from AMY's integral equations as in Eq. (\ref{eq:1to2}). 

\subsubsection{Treatment of hard interactions: elastic case ($2\leftrightarrow 2$)}

\label{sec:factorized_model:elastic_hard}

\begin{figure}
    \centering
    \includegraphics[width=0.45\linewidth]{plots/large_angle_diagram.pdf}
    \includegraphics[width=0.45\linewidth]{plots/splitting_diagram.pdf}
    \caption{(a) Elastic $2\leftrightarrow 2$ interactions with $\tilde{q}_\perp > \mu_{\tilde{q}_\perp}$ and $\omega < \Lambda$ is denoted as large-angle interaction; (b) Elastic $2\leftrightarrow 2$ interactions with $\omega > \Lambda$ are treated with a splitting approximation (see text). }
    \label{fig:elas_diagram}
\end{figure}
The diagram of a typical $2\leftrightarrow 2$ elastic interaction is shown in Fig. \ref{fig:elas_diagram}. We define the momentum transfer between the two incoming particles as $Q=(\omega, \vec{q})$; the four-momenta of the incoming and outgoing particles are $P=(p, \mathbf{p})$, $K=(k, \mathbf{p})$, $P' = P-Q$, and $K'=K+Q$. Using $\tilde{q}_\perp \equiv \sqrt{q^2-\omega^2}$ and $\omega$,  we divide the phase space of elastic interactions as
\begin{itemize}
    \item Large-angle scattering $\mathcal{C}_{\textrm{large-angle}}^{2\leftrightarrow2}$:  $\tilde{q}_\perp > \mu_{\tilde{q}_\perp}$ and $\omega < \Lambda$;
    \item ``Splitting-like'' process  $\mathcal{C}_{\textrm{split}}^{2\leftrightarrow2}$: $\omega > \Lambda$
\end{itemize}
with
\begin{equation}
\mathcal{C}^{2\leftrightarrow 2}_{\textrm{hard}} = \mathcal{C}_{\textrm{large-angle}}^{2\leftrightarrow2}+\mathcal{C}_{\textrm{split}}^{2\leftrightarrow2}.
\end{equation}

\paragraph{Large-angle scattering ($\mathcal{C}_{\textrm{large-angle}}^{2\leftrightarrow2}$)}

Hard scatterings with $\tilde{q}_\perp > \mu_{\tilde{q}_\perp}$ and $\omega < \Lambda$ are denoted as large-angle interactions, because the scattering angle 
\begin{equation}
\frac{q_\perp}{q_z} \approx \frac{\tilde{q}_\perp}{\omega} 
\end{equation}
is generally large in this region. The cutoff $\mu_{\tilde{q}_\perp}$ is typically assumed to be $gT \ll \mu_{\tilde{q}_\perp} \ll T$~\cite{Ghiglieri:2015ala}, although we will see in Section~\ref{sec:tequila_weak_coupling} that this condition can be relaxed. We assume $p \gg \omega$ in this region, and simplify the matrix elements accordingly. 

We use vacuum matrix elements for $\mathcal{C}^{2\leftrightarrow 2}_{\textrm{hard}}$, because the screening effects are only significant for soft interactions ($\mathcal{C}^{2\leftrightarrow 2}_{\textrm{soft}}$) in the weakly-coupled regime. Since we are only interested in the evolution of energetic partons, we keep terms to the first order in $T/p$ in the matrix elements. 

The treatment of the region $\tilde{q}_\perp > \mu_{\tilde{q}_\perp}$ and $p-\Lambda < \omega < p$ --- which is handled differently for technical reasons --- is discussed in \app{app:remaining}. 

\paragraph{Splitting approximation ($\mathcal{C}_{\textrm{split}}^{2\leftrightarrow2}$)}
When both of the outgoing particles of a $2\leftrightarrow 2$ interaction are hard ($p', k' > p_{\textrm{cut}}$), the interaction can be effectively considered as a splitting process. The splitting leads to a hard recoil parton which should be included in the calculation. 

We use a cutoff $\Lambda$ on $\omega$ to distinguish two hard outgoing particles from only one hard outgoing particles. In principle, this cutoff $\Lambda$ should be $3T \ll \Lambda \ll p$. In the numerical implementation, unless specified otherwise, we choose $\Lambda = \min(\sqrt{3pT}, p_{\textrm{cut}})$ to divide $\mathcal{C}^{2\leftrightarrow 2}_{\mathrm{split}}$ and $\mathcal{C}^{2\leftrightarrow 2}_{\mathrm{large-angle}}$. Recall that we use $p_{\textrm{cut}}=2$~GeV in this work. As shown in Fig. \ref{fig:elas_diagram}, splitting approximation process is the $2\leftrightarrow 2$ interactions with $\Lambda < \omega < p-\Lambda$.  

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{plots/hard_elas_differential_rate.pdf}
    \caption{The differential rate of splitting approximation interactions and large-angle interactions for $gg \leftrightarrow gg$ process when $\alpha_s=0.3$. The shaded area is the region of $\sqrt{3p_{\textrm{cut}}T} < \omega < p_{\textrm{cut}}$. $d\Gamma_{\textrm{vac}}/d\omega$ is the differential rate of vacuum matrix elements for $2\leftrightarrow 2$ interactions. The results are for $p_0=100$~GeV. Note that in the numerical implementation, we double-count the large-angle interaction rate because we only sample in half of the phase space. Here, to compare with splitting approximation rate, we decrease the large-angle interaction rate in the numerical implementation by a factor of $\frac{1}{2}$ to cancel out the double-count. }
    \label{fig:hard_elas}
\end{figure}

At the interface between the phase space of large-angle scattering ($\mathcal{C}_{\textrm{large-angle}}^{2\leftrightarrow2}$) and splitting-like processes ($\mathcal{C}_{\textrm{split}}^{2\leftrightarrow2}$), the two collision kernels should be consistent.
We verified this in Fig.~\ref{fig:hard_elas}: the differential rates of $\mathcal{C}_{\textrm{split}}^{2\leftrightarrow2}$ and $\mathcal{C}_{\textrm{large-angle}}^{2\leftrightarrow2}$ are compatible in $\sqrt{3p_{\textrm{cut}} T} < \omega < p_{\textrm{cut}}$. As long as we choose the cutoff $\Lambda$ in this range, this division of $\mathcal{C}^{2\leftrightarrow 2}_{\textrm{hard}}$ should be consistent.

A detailed discussion of the splitting approximation process is in \app{app:splittingapprox}. 
The  $p \gg T$ and $\omega \gg T$ kinematic cuts lead to significant simplifications for the matrix elements entering into $\mathcal{C}_{\textrm{split}}^{2\leftrightarrow2}$.



\subsubsection{Treatment of soft interactions}

\label{sec:hard_soft_theory:reformulation:soft}


In the hard-soft factorized approach, the large number of  soft interactions are described stochastically with drag and diffusion coefficients. When the momentum transfer is small, the Boltzmann equation [Eq. (\ref{eq:Boltzmann})] can be approximated as a Fokker-Planck equation. The collision kernel of the Fokker-Planck equation is written as: 
\begin{eqnarray}
\mathcal{C}_{\textrm{diff}}^{1\leftrightarrow 2, 2\leftrightarrow2} & = & \mathcal{C}^{1\leftrightarrow 2}_{\textrm{soft}}[\delta f]+\mathcal{C}^{2\leftrightarrow 2}_{\textrm{soft}}[\delta f] \nonumber \\
&= & -\frac{\partial}{\partial p^i}\left[\eta_{D, \textrm{soft}}p^i\delta f\right] \nonumber \\
& & -\frac{1}{2}\frac{\partial^2}{\partial p^i\partial p^j}\left\{\left[\hat{p}^i\hat{p}^j \hat{q}_L+\frac{1}{2}\left(\delta^{ij}-\hat{p}^i\hat{p}^j\right)\hat{q}\right]\delta f\right\} \nonumber \, , \\
    \label{eq:FP}
\end{eqnarray}
where $\eta_{D, \textrm{soft}}$ is the drag coefficient of the soft interactions, $\hat{q}_{L, \textrm{soft}}$ and $\hat{q}_\textrm{soft}$ are the longitudinal and transverse momentum diffusion coefficients of the soft interactions. 

In the diffusion process, the number and the identity of the particles are preserved. Since the soft radiations of the $1 \leftrightarrow 2$ process are absorbed by the plasma, the number of particles is also preserved.
We include both the soft $1\leftrightarrow 2$ and $2\leftrightarrow 2$ collisions in the diffusion process. The diffusion process can be solved using a Langevin equation~\cite{He:2013zua} in the numerical implementation.

For soft $1\leftrightarrow 2$ process, we can obtain the perturbative longitudinal diffusion coefficients by expanding $\mathcal{C}^{1\leftrightarrow 2}$ and only keeping the soft radiation terms. At  leading order in $\alpha_s$,
\begin{equation}
\hat{q}^{1\leftrightarrow 2}_{L, \textrm{soft}} = \frac{(2-\ln2)g^4C_RC_AT^2\mu_\omega}{4\pi^3} \, ,
\label{eq:inel_qhat_L}
\end{equation}
where $C_R$ is the Casimir factor.\footnote{Note that the diffusion coefficient $\hat{q}^{1\leftrightarrow 2}_{L, \textrm{soft}}$ does not depend on the number of the quark flavor, because very soft radiations are dominated by gluon scatterings.} The derivation of this value can be found in Appendix~\ref{appendix:inel_low_omega}. 

We assume that the radiation angle is zero for collinear radiations. Consequently the transverse diffusion coefficient of $1\leftrightarrow 2$ interactions is approximated as zero. 


For soft $2\leftrightarrow 2$ processes, 
the diffusion coefficients can be calculated perturbatively;
a modern derivation  
can be found in Ref.~\cite{Ghiglieri:2015ala}.
The transverse momentum diffusion coefficient due to soft scatterings is
\begin{equation}
   \hat{q}_{\rm soft}^{2 \leftrightarrow 2} 
   = \frac{g^2 C_R T m_D^2}{4\pi} \ln \left[1+\left(\frac{\mu_{\tilde{q}_\perp}} {m_D}\right)^2\right] \, ,
    \label{eq:elas_qhat}
 \end{equation}
where 
$m_D^2 \equiv g^2T^2(N_c/3+N_f/6)$ is the square of the leading order Debye mass, $N_c=3$ is the number of colors and $N_f$ is the number of flavors involved in the interactions. 
The longitudinal diffusion coefficient at order $\mathcal{O}(\alpha_s)$ is
\begin{equation}
    \hat{q}^{2\leftrightarrow 2}_{L, \textrm{soft}} = \frac{g^2C_RTM_\infty^2}{4\pi}\ln\left[1+\left(\frac{\mu_{\tilde{q}_\perp}}{M_\infty}\right)^2\right]\, ,
    \label{eq:elas_qhat_L}
\end{equation}
where $M_\infty \equiv \sqrt{m_D^2/2}$ is the gluon asymptotic thermal mass. 
For both diffusion coefficients, the cutoff $\mu_{\tilde{q}_\perp}$ is assumed to be of order $T$ or smaller.

Since detailed balance is preserved in the Fokker-Planck equation, as verified in  \app{app:detailed_balance}, the drag coefficient $\eta_D$ can be calculated from diffusion coefficients according to Einstein relation for both soft $1\leftrightarrow 2$ and $2\leftrightarrow 2$ processes: 
\begin{equation}
\eta_{D, \textrm{soft}}(E) = \frac{\hat{q}_{L, \textrm{soft}}}{2Tp}\left[1+\mathcal{O}\left(\frac{T}{p}\right)\right]\, . 
\label{eq:drag}
\end{equation}

Equations (\ref{eq:inel_qhat_L}-\ref{eq:drag}) are the perturbative drag and diffusion coefficients derived in the weakly coupled regime. Note that these perturbative transport coefficients are not necessarily valid in the phenomenological regime studied in heavy ion collisions. We discuss the range of validity of the perturbative coefficients in Section \ref{sec:weak_coupling_analytical_vs_numerical_qhat}. 

Besides the identity preserving diffusion process, the identity of the particle can be converted through soft fermion
exchange with the medium. This exchange must be screened with the non-perturbative HTL resummation scheme. In the hard-soft factorized approach adopted here, we separate the $2\leftrightarrow 2$ processes with fermion exchange into 
hard  collisions with $\tilde q_\perp > \mu_{\tilde q_\perp}$,  and soft collisions with $\tilde q_\perp < \mu_{\tilde q_\perp}$ (see \Fig{fig:reformulation}). The hard exchange collisions are treated with vacuum matrix elements, while the soft exchange collisions are incorporated into  
a conversion rate $\Gamma_{q \rightarrow g}^{\rm conv}(p)$ for $q \rightarrow g$:
\begin{equation}
\Gamma_{q\rightarrow g}^{\rm conv}(p) =  \frac{g^2C_F m_\infty^2 }{16\pi p} \log\left[1 + \frac{\mu_{\tilde q_\perp}^2}{m_{\infty}^2} \right] \, .
\end{equation}
Here $m_{\infty}^2$ is the fermion asymptotic mass, $m_{\infty}^2 = g^2 C_F T^2/4$.
In each time step there is a probability $\Delta t \, \Gamma^{\rm conv}$ for a quark to become a gluon, with the same momentum, and vice versa.
Further details about the conversion rate $C_{\rm conv}^{2\leftrightarrow 2}$ are given in \app{app:soft_conv}.
In the future, the non-perturbative conversion coefficient $\Gamma_{q \rightarrow g}^{\rm conv}$ can be taken from 
a next-to-leading order analysis~\cite{Ghiglieri:2015ala}, 
or can be determined from data in a Bayesian approach.

\subsubsection{Summary}

In summary, the collision kernel of hard-soft factorized model is reformulated as
\begin{eqnarray}
\mathcal{C} &=& \mathcal{C}^{2\leftrightarrow 2} + \mathcal{C}^{1\leftrightarrow 2} \nonumber \\
&=& \mathcal{C}_{\textrm{large-}\omega}^{1\leftrightarrow2}\left(\mu_\omega\right)+\mathcal{C}_{\textrm{large-angle}}^{2\leftrightarrow2}\left(\mu_{\tilde{q}_\perp}, \Lambda\right)+ \mathcal{C}_{\textrm{split}}^{2\leftrightarrow2}\left(\Lambda\right) \nonumber \\
& &\;\;\;\; +\mathcal{C}_{\textrm{diff}}^{1\leftrightarrow2, 2\leftrightarrow2}\left(\mu_\omega, \mu_{\tilde{q}_\perp}\right)+\mathcal{C}^{2\leftrightarrow 2}_{\textrm{conv}}\left(\mu_{\tilde{q}_\perp}\right) \, .
\label{reformulation}
\end{eqnarray}

The cutoff dependence of the stochastic description is cancelled in Eq.~(\ref{reformulation}) by the cutoff dependence of the hard interactions.
That is, each individual process in the hard-soft factorized model is dependent on the cutoff, but this dependence cancels out when all the processes are summed. 
We show this explicitly in Section~\ref{sec:tequila_weak_coupling}. 



\subsection{Running of the strong coupling $\alpha_s$}

All discussions up to this point assumed that the strong coupling constant $\alpha_s$ is fixed at a given small value. 
It is clear, however, that the strong coupling constant will be different for soft and hard interactions; this is in fact a key assumption of the present model: hard interactions are more perturbative than soft ones, because the coupling constant scales inversely with the momentum exchange between the energetic parton and the plasma (see Ref.~\cite[Section V]{CaronHuot:2010bp} for a discussion, for example).
The running is slow (logarithmic in the momentum exchange), however, more studies will be necessary to understand the exact magnitude of loop corrections or non-perturbative effects on soft and hard collisions.

As a first step in introducing our model of parton energy loss, we keep the strong coupling constant $\alpha_s$ fixed throughout the manuscript.

\section{Hard-soft factorization of parton energy loss in the weakly-coupled regime: numerical study}

\label{sec:tequila_weak_coupling}

In the first part of this section, we compare the analytical equations for the soft-interaction parton transport coefficients [Eqs.~(\ref{eq:inel_qhat_L}--\ref{eq:elas_qhat_L})] with their numerical values evaluated from the matrix elements, and summarize the range of cutoff and coupling where they are consistent. We also compare (i) soft interactions modelled with matrix elements with (ii) soft interactions modelled with the Langevin equation. We perform this test in the weak coupling limit. We use this discussion to review the range of validity of the Fokker-Planck equation and its stochastic Langevin implementation. 

In the second part of this section, we compute the energy loss of an energetic parton in a brick and discuss the dependence of the results on the soft-hard cutoffs introduced in Section~\ref{sec:factorized_model}.


\subsection{Treatment of soft interactions}

Soft interactions can be described either stochastically with transport coefficients, or microscopically with matrix elements.
In what follows, we compare these two descriptions,
with particular emphasis on the effect of the soft-hard cutoffs and of the coupling constant.
All calculations in this subsection are performed in the pure glue limit ($N_f=0$) for simplicity. 



\subsubsection{Analytical and numerical calculation of soft transport coefficients}
\label{sec:weak_coupling_analytical_vs_numerical_qhat}


In a weakly-coupled quark-gluon plasma, the  drag and diffusion coefficients for soft interactions can be calculated analytically using perturbation theory [~Eqs.(\ref{eq:inel_qhat_L}-\ref{eq:elas_qhat_L})~], as discussed in Section~\ref{sec:hard_soft_theory:reformulation:soft}.
The same drag and diffusion coefficients can be obtained by direct numerical integration of the parton energy loss rates; these rates are calculated from matrix elements screened by plasma effects~\cite{Arnold:2003zc}.

The diffusion coefficients are defined as \cite{Ghiglieri:2015ala}
\begin{equation}
\begin{split}
    \hat{q}(p) &\equiv \frac{d}{dt}\left<\left(\Delta p_\perp\right)^2\right>, \\
    \hat{q}_L(p) &\equiv \frac{d}{dt}\left<\left(\Delta p_L\right)^2\right>,
\end{split}
    \label{eq:df_qhat}
\end{equation}
where $\Delta p_\perp$ is the momentum change perpendicular to the direction of the energetic parton, and $\Delta p_L$ is the longitudinal momentum change of the parton. The brackets represent an average over all interactions during the parton propagation. The numerical soft diffusion rates are thus calculated as 
\begin{equation}
\begin{split}
    \hat{q}^{2\leftrightarrow 2}_{\textrm{soft}}(p) &= \int_0^{\mu_{\tilde{q}_\perp}} d\tilde{q}_\perp \int_{-\infty}^{\Lambda} d\omega \tilde{q}_\perp^2 \left.\frac{d^2\Gamma(\bm{p}, \bm{q})}{d\omega d\tilde{q}_\perp}\right|_{2 \leftrightarrow 2}, \\
    \hat{q}^{2\leftrightarrow 2}_{L, \textrm{soft}}(p) &= \int_0^{\mu_{\tilde{q}_\perp}}d\tilde{q}_\perp \int_{-\infty}^{\Lambda} d\omega \omega^2\left.\frac{d^2\Gamma(\bm{p}, \bm{q})}{d\omega d\tilde{q}_\perp}\right|_{2 \leftrightarrow 2}, \\
    \hat{q}^{1\leftrightarrow 2}_{L, \textrm{soft}}(p) &= \int_{-\mu_\omega}^{\mu_\omega} d\omega \omega^2\left.\frac{d\Gamma(\bm{p}, \bm{q})}{d\omega}\right|_{1 \leftrightarrow 2},
\end{split}
\label{eq:num_qhat}
\end{equation}
where $d\Gamma(\bm{p}, \bm{q})/d\omega$ and $d^2\Gamma(\bm{p}, \bm{q})/d\omega d\tilde{q}_\perp$ are the rates for an energetic parton with four-momentum $(p,\bm{p})$ to undergo a four-momentum change $(\omega,\bm{q})$ calculated using screened matrix elements. The initial parton energy $p$ is assumed to be much larger than all other energy scales in the problem, effectively $p\to\infty$.
The cutoffs $\mu_{\tilde{q}_\perp}$, $\Lambda$ and $\mu_\omega$ are used to limit the phase space of interactions included in the transport coefficients, in the present case to limit the interactions to soft ones only.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{plots/qhat_ratio_cutoff_largerange.pdf}
    \caption{The ratio between the numerical (\eqref{eq:num_qhat}) and analytical (Eqs.(\ref{eq:inel_qhat_L}-\ref{eq:elas_qhat_L})) 
    momentum transport coefficients: $\hat{q}^{1\leftrightarrow 2}_{L, \textrm{soft}}$, $\hat{q}^{2\leftrightarrow 2}_{\textrm{soft}}$ and $\hat{q}^{2\leftrightarrow 2}_{L, \textrm{soft}}$. The numerical
    results are computed 
   
    with exact $1\leftrightarrow 2$ or $2\leftrightarrow 2$ kinematics up to a cutoff $\mu/T$. The analytical coefficients make kinematic approximations appropriate for $\mu/T \ll 1$.
    The results are shown for different values of the hard-soft cutoffs at $\alpha_s=0.005$. We calculate these results using $p_0 = 100$~GeV and $T=300$~MeV. The cutoff $\mu$ in the figure denotes $\mu_{\tilde{q}_\perp}$ for $\hat{q}^{2\leftrightarrow 2}_{\textrm{soft}}$ and $\hat{q}^{2\leftrightarrow 2}_{L, \textrm{soft}}$, and $\mu_\omega$ for $\hat{q}^{1\leftrightarrow 2}_{L, \textrm{soft}}$. In the elastic case, the additional cutoff on $\omega$ is set to $\Lambda=\min(p_{\textrm{cut}}, \sqrt{3p_0T})$. }
    \label{fig:qhat_cutoff}
\end{figure}

There are two important differences between Eq.~(\ref{eq:num_qhat}) and the analytical diffusion coefficients Eqs.~(\ref{eq:inel_qhat_L}-\ref{eq:elas_qhat_L}). 
First, Eq.~(\ref{eq:num_qhat}) is formally valid for arbitrarily large cutoffs ($\mu_{\tilde{q}_\perp}$, $\Lambda$ and $\mu_{\omega}$), while Eqs.~(\ref{eq:inel_qhat_L}- \ref{eq:elas_qhat_L}) assume the cutoff to be at most of order $T$.
Second, there is the question of the smallness of the coupling. Equations~(\ref{eq:inel_qhat_L}- \ref{eq:elas_qhat_L}) are derived assuming $\alpha_s\ll 1$. Equation~(\ref{eq:num_qhat}) is valid at arbitrarily coupling, although the rates $d\Gamma(\bm{p}, \bm{q})/d\omega$ and $d^2\Gamma(\bm{p}, \bm{q})/d\omega d\tilde{q}_\perp$ themselves are typically calculated perturbatively.\footnote{It is highlighted in Ref.~\cite{CaronHuot:2010bp} that the AMY differential equation used to evaluate the inelastic collisions rate remains similar if interactions with the plasma are non-perturbative. One difference is the perturbative partonic collision kernel $C(\mathbf{q}) \propto m_D^2 /\left[ \mathbf{q}^2 (\mathbf{q}^2+m_D^2) \right]$ that must be modified. Non-perturbative contributions to the thermal masses are another difference.\label{footnote:caronhuot}}

A comparison of Eq.~(\ref{eq:num_qhat}) and the analytical diffusion coefficients Eqs.~(\ref{eq:inel_qhat_L}-\ref{eq:elas_qhat_L}) is shown in Fig.~\ref{fig:qhat_cutoff} as a function of the different cutoffs. This comparison is made at weak coupling ($\alpha_s=0.005$) and yields the expected agreement between the two approaches, as long as the cutoffs are $\lesssim T$.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{plots/qhat_compare_coupling.pdf}
    \caption{Comparison of the numerical and analytical $\hat{q}^{1\leftrightarrow 2}_{L, \textrm{soft}}$, $\hat{q}^{2\leftrightarrow 2}_{\textrm{soft}}$ and $\hat{q}^{2\leftrightarrow 2}_{L, \textrm{soft}}$ with different coupling constants $\alpha_s$ (see Fig.~\ref{fig:qhat_cutoff} for description). For the kinematic cutoffs, we use $\mu_{\tilde{q}_\perp} = \mu_\omega = T$ and $\Lambda=\min(\sqrt{3p_0T}, p_{\textrm{cut}})$. The numerical values of the transport coefficients were calculated assuming $T=300$~MeV and an energetic parton with $p_0 = 100$~GeV. }
    \label{fig:qhat_coupling}
\end{figure}

In Fig.~\ref{fig:qhat_coupling} we compare the analytical soft diffusion coefficients Eqs.~(\ref{eq:inel_qhat_L}-\ref{eq:elas_qhat_L}) with the numerical soft diffusion coefficients Eq.~(\ref{eq:num_qhat}) at different values of the strong coupling constant $\alpha_s$. 
We find that the analytical soft diffusion coefficients agree well with the numerical calculations even at large coupling, except for a small tension in $\hat{q}^{2\leftrightarrow 2}_{L, \textrm{soft}}$ at large $\alpha_s$.
Tension between different calculations of the soft transport coefficients are in fact not unexpected: perturbative calculations can be equivalent at order $g^n$ yet be different at order $g^{n+1}$. These differences are negligible at weak coupling, but can become significant for larger values of the coupling. This is a natural consequence of pushing the calculations beyond their regimes of validity. 
There is a practical consequence: two parton energy loss calculations that use the exact same approach (weakly-coupled kinetic theory) can lead to different results, when used at large coupling; neither approach is more ``correct'' than the other. This is important to keep in mind when comparing the present soft-hard factorized energy loss model with other implementations such as Ref.~\cite{Schenke:2009gb}.



\subsubsection{Theoretical guidance on the range of applicability of the Fokker-Planck equation}

\label{sec:hard_soft_at_weak_coupling:soft_interactions:scale}

The energy loss of energetic partons through soft interactions is described by solving the Fokker-Planck equation with a stochastic Langevin approach. The applicability of the stochastic description is limited to the regime where the Fokker-Planck equation holds. 
This regime of applicability depends partly on properties of the interactions rates.
We can summarize the regime of validity of the Fokker-Planck equation by first expanding the \emph{Boltzmann equation} for soft collisions (around $\omega=0$): 
\begin{equation}
\begin{split}
    \partial_t f(p, t) = &\left<\omega\right> f^{(1, 0)}(p, t) +\frac{1}{2}\left<\omega^2\right>f^{(2, 0)}(p, t)\\ 
    &+\frac{1}{6}\left<\omega^3\right>f^{(3, 0)}(p, t)+\dots, 
\end{split}
\label{eq:boltzmann_small_w_expansion}
\end{equation}
where $f(p,t)$ is the momentum distribution of energetic partons at time $t$ and
\begin{equation}
    \left<\omega^k\right> =  \int d\omega \omega^k\frac{d\Gamma}{d\omega}
    \label{eq:rate_moment}
\end{equation}
is the k-th moment of the differential collision rate $d\Gamma/d\omega$.\footnote{
The bounds on the integration are the same as in Eqs.~(\ref{eq:num_qhat}), including the additional integration over $\tilde{q}_\perp$ necessary in the elastic case.}


By keeping only the first two terms on the right-hand side, \eqref{eq:boltzmann_small_w_expansion} simplifies to the Fokker-Planck equation. 

Assuming a single initial energetic parton of energy $p_0$,
\begin{equation}
    f(p,t=0)=\delta(p-p_0),
\end{equation}
the solution of the Fokker-Planck equation is
\begin{equation}
    \begin{split}
        f_{FP}(p, t)&=\frac{\exp{\left[-\frac{\left( p-(p_0-\left<\omega\right>t ) \right)^2 }{2t\left<\omega^2\right>}\right]}}{\sqrt{2\pi\left<\omega^2\right> t}} \, .
    \end{split}
    \label{eq:toy_FP}
\end{equation}
The above solution simply describes the energy distribution of the energetic parton widening from scatterings with $\hat{q}_L = \left<\omega^2\right>$ energy diffusion, and an average energy loss of $\left<\omega\right> t$.

Using this solution, we can compute the ratio of the third and second terms in the expanded Boltzmann equation (\eqref{eq:boltzmann_small_w_expansion}): 
\begin{equation}
    \mathcal{R} = \frac{\frac{1}{6}\left<\omega^3\right>f_{FP}^{(3, 0)}(p, t)}{\frac{1}{2}\left<\omega^2\right>f_{FP}^{(2, 0)}(p, t)} = -\frac{2\Delta p \left<\omega^3\right>(\Delta p^2-3\left<\omega^2\right>t)}{3\left<\omega^2\right>t(\Delta p^2-\left<\omega^2\right>t)} \, ,
    \label{eq:R_factor}
\end{equation}
where $\Delta p = p-(p_0-\left<\omega\right>t)$ is the distance in momentum from the peak of the Fokker-Planck solution (\eqref{eq:toy_FP}). 

Significant corrections to the Fokker-Planck solution \eqref{eq:toy_FP} are expected unless $\mathcal{R} \ll 1$.
As is clear from \eqref{eq:R_factor}, the range of validity of the Fokker-Planck equation depends on properties of the rate (the second and third moments $\left<\omega^2\right>$ and $\left<\omega^3\right>$), as well as on time $t$ and on the distance in momentum $\Delta p$ from the peak of the distribution.

The Fokker-Planck equation describes the effect of soft interactions on an energetic parton. 
The soft interactions dominate for small values of $\Delta p$.
Expanding $\mathcal{R}$ (\eqref{eq:R_factor}) around $\Delta p=0$, we obtain:
\begin{equation}
    \mathcal{R} = -\frac{\Delta p \left<\omega^3\right>}{\left<\omega^2\right>^2 t}+\frac{2\Delta p^3\left<\omega^3\right>}{3\left<\omega^2\right>^3t^2}+\mathcal{O}(\Delta p^5) \, .
    \label{eq:R_factor_expansion}
\end{equation}

By taking the ratio of the second and first term of this expansion, 
\begin{equation}
r \equiv \frac{-2\Delta p^3\left<\omega^3\right>}{3\left<\omega^2\right>^3t^2} \bigg/ \frac{\Delta p \left<\omega^3\right>}{\left<\omega^2\right>^2 t},
\end{equation}
we can find the value of $\Delta p$ for which this ratio will be large:
\begin{equation}
\Delta p = \sqrt{\frac{3}{2}r}\sqrt{\left<\omega^2\right>t} \, ,
\label{eq:Delta_p}
\end{equation}
with $r$ a constant assumed to be smaller than $1$.
We can use this value of $\Delta p$ as the range of momentum around the mean energy loss that can reasonably be described by the Fokker-Planck equation. 
Using \eqref{eq:Delta_p} and the first term of \eqref{eq:R_factor_expansion}, we define the scale $\mathcal{S}$ as
\begin{equation}
    \mathcal{S}=\frac{\left<\omega^3\right>}{\left<\omega^2\right>^{3/2}} \frac{1}{\sqrt{t}} \, .
    \label{eq:scale}
\end{equation}

When this scale $\mathcal{S}$ is much smaller than 1, the Fokker-Planck equation is expected to provide a good description of the Boltzmann equation in the relevant range of momentum. 

For inelastic interactions at weak coupling, the scale is given by
\begin{equation}
\mathcal{S}_{1\leftrightarrow 2}=\frac{\pi ^{3/2}}{3 C_A \sqrt{2-\ln(2)}} \frac{\mu^{3/2}}{g^2  T^2 \sqrt{t}}.
\label{eq:scale_inel}
\end{equation}
This implies that soft inelastic emissions with energy smaller than $\mu$ can be described with the Langevin equation as long as the evolution time $t$ in the medium is sufficiently long:
\begin{equation}
    t \gg \frac{\mu^3}{g^4 T^4}.
    \label{eq:time_langevin_inel}
\end{equation}
Assuming $\mu \sim T$ results in $t \gg 1/[g^4 T]$, while $\mu \sim g T$ results in $t \gg 1/[g T]$. 
This implies that there is a very large difference between a stochastic description of soft interactions with $\omega \lesssim T$ compared to soft interactions with $\omega \lesssim g T$: in the former case, one needs a plasma $1/g^3$ larger.
These values serve as a reminder that, while one can in principle increase the phase space of interactions described stochastically,
one may need an unrealistically large plasma for this description to be valid.




\subsubsection{Comparison between the diffusion process and the collision rate}

\label{sec:weak_coupling_diffusion_vs_rate}

\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{plots/soft_collRate_and_langevin_comparison_smallcoupling_1fm.pdf}
    \caption{The energy distribution of a $100$~GeV gluon propagating through a $300$~MeV pure glue medium at $\alpha_s = 0.005$. The evolution time is $t = (0.3/\alpha_s )^2=3600$~fm. Only soft $1\leftrightarrow 2$ interactions with $\omega < \mu_\omega$ are allowed. Three different values of the cutoff are shown: $\mu_\omega/T=0.5, 1, 2$.}
    \label{fig:diffusion_collision_rate_small_coupling:inel}
\end{figure*}

\begin{figure*}
    \includegraphics[width=2\columnwidth]{plots/soft_collRate_and_langevin_comparison_elas_small_coupling_numerical_1fm_scale.pdf}
    \onecolumngrid
    \caption{The energy distribution of a $100$~GeV gluon propagating through a $300$~MeV pure glue medium at $\alpha_s = 0.005$. The evolution time is $t = (0.3/\alpha_s )^2=3600$~fm. Only soft $2\leftrightarrow 2$ interactions with $\omega < \Lambda$ and $\tilde{q}_\perp < \mu_{\tilde{q}_\perp}$ are allowed. Three different values of the cutoff $\mu_{\tilde{q}_\perp}$ are shown: $\mu_{\tilde{q}_\perp}/T=0.5, 1, 2$. }
    \label{fig:diffusion_collision_rate_small_coupling:elastic}
\end{figure*}

In this section, we verify numerically the conclusion from the previous section:
we compare a stochastic and a microscopic evolution of energetic partons in a static medium. In the microscopic rate-based picture, we use kinematic cuts to forbid hard interactions of the energetic parton. 
Because we are comparing soft interactions, we must use screened elastic matrix elements~\cite{Arnold:2003zc} in the microscopic description.\footnote{Note that this is for testing purpose only, and that this is different from the vacuum matrix elements used for $\mathcal{C}^{2\leftrightarrow 2}_{\textrm{hard}}$ in the hard-soft factorized energy loss model.} 
The screened inelastic ($1\leftrightarrow 2$) rate is obtained numerically by solving the AMY differential equation, except for very small $\omega$ values, in which case the analytical expression described in Appendix~\ref{appendix:inel_low_omega} (\eqref{eq:inel_soft_rate}) is used.

We choose the hard-soft cutoffs (i.e. $\mu_\omega$ and $\mu_{\tilde{q}_\perp}$) to be at the order of $T$ in the following tests.
We set the coupling to be $\alpha_s = 0.005$, which corresponds to $g\approx 0.25$.
We choose $T=300$~MeV for the temperature of the fluid, and set the propagation time in the plasma to be $t=(0.3/\alpha_s)^2=3600$~fm.\footnote{We choose the evolution time $t \propto 1/\alpha_s^2$ to keep the number of the collisions approximately the same for different values of $\alpha_s$. With the choice $t=(0.3/\alpha_s)^2$, the evolution time is $1$~fm when we use $\alpha_s=0.3$ later in the manuscript.}


We perform the diffusion approach and the collision rate approach separately to calculate the single parton energy distribution of a hard $100$~GeV gluon propagating in the static pure glue medium. We emphasize once again that we only include soft interactions in the test by introducing the following hard-soft cutoffs on radiation energy and momentum transfer: for $\mathcal{C}^{1\leftrightarrow 2}_{\textrm{soft}}$, we only include radiations with the radiation energy $\omega < \mu_\omega$; while for $\mathcal{C}^{2\leftrightarrow 2}_{\textrm{soft}}$, we only include interactions with $\tilde{q}_\perp < \mu_{\tilde{q}_\perp}$ and the energy transfer $\omega < \Lambda$. 

According to \eqref{eq:scale}, for inelastic interactions ($\mathcal{C}^{1\leftrightarrow 2}_{\textrm{soft}}$) to be describable stochastically for a cutoff $\sim T$, one needs $t \gg 1/[g^4 T] \approx 200$~fm of propagation time in the conditions described above. As expected, we find in Figure \ref{fig:diffusion_collision_rate_small_coupling:inel} that for inelastic interactions, in the weakly-coupled regime, the diffusion process can reproduce the single parton energy distribution generated by the collision-rate process. The value of the scale $\mathcal{S}$, shown for each cutoff $\mu_\omega$, are indeed smaller than $1$. As $\mu_\omega$ increases, small differences appear between the Langevin description and the microscopic collision approach; the scale $\mathcal{S}$ is correspondingly larger, though still smaller than $1$.

The same results are shown for the elastic case ($\mathcal{C}^{2\leftrightarrow 2}_{\textrm{soft}}$) in Fig.~\ref{fig:diffusion_collision_rate_small_coupling:inel}.  This time, the scale $\mathcal{S}$ is somewhat larger, and somewhat larger differences can indeed be seen between the Langevin and collision rate descriptions. As for the elastic case, the scale $\mathcal{S}$ increases as the cutoff increases, where more and more collisions are described stochastically.


\subsection{Parton energy loss at small coupling in a static medium}

\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{plots/cutoff_dependence_small_coupling.pdf

    \caption{The energy distribution of a $100~\textrm{GeV}$ gluon propagating through $300~\textrm{MeV}$ QGP medium ($N_f=3$) at $\alpha_s = 0.005$ with different values of the cutoffs. The evolution time is $t=(0.3/\alpha_s )^2=3600$~fm. The subplot (a) only includes the $\mathcal{C}^{1\leftrightarrow 2}$ interactions and (b) only includes $\mathcal{C}^{2\leftrightarrow 2}$ interactions. In both cases the cutoff $\mu$ is varied: the soft interactions (those with momentum transfer less than $\mu_\omega$ and $\mu_{\tilde{q}_\perp}$ respectively) are treated with a Langevin process, while the rest of the kinematic phase space is treated with rates. }
    \label{fig:cutoff_dependence_small_coupling}
\end{figure*}

Building on the validation from the previous section,
we can combine our approaches for the hard and soft interactions to implement the entire hard-soft factorized parton energy loss model described in Section~\ref{sec:factorized_model}. Remember that in the following, we use vacuum matrix elements for $\mathcal{C}^{2\leftrightarrow 2}_{\textrm{hard}}$, since the screening effects are encoded in the drag and diffusion coefficients of soft interactions. We also extend this test to a full quark-gluon plasma, with $N_f=3$. We use once again $\alpha_s = 0.005$ ($g\approx 0.25$), with a propagation time of $t=(0.3/\alpha_s)^2=3600$~fm in a $T=300$~MeV plasma.


As summarized by Eq. (\ref{reformulation}), the hard or soft processes alone are dependent on the cutoff, but their cutoff dependence cancels out when combined. 
We confirm that, for both the inelastic and elastic cases, the single parton energy distribution is independent on the hard-soft cutoffs at small coupling
in Fig. \ref{fig:cutoff_dependence_small_coupling}, given a sufficiently long evolution time. These results are consistent with those obtained in the previous section.

\section{Hard-soft factorization of parton energy loss beyond weak-coupling}
\label{sec:large_coupling}

Soft interactions between an energetic parton and a deconfined plasma are likely non-perturbative. Evaluating this non-perturbative rate from first principles is an ongoing challenge. In this section, we estimate this non-perturbative rate using a typical approach in the heavy ion literature: we use the perturbative rate and extrapolate it to large coupling. 

Recall that we do not use a running coupling in this work. As such, we use the same value of $\alpha_s$ for soft and hard interactions, with the understanding that the future introduction of a running coupling will indeed lead to smaller values of $\alpha_s$ for hard interactions, as assumed in this work.

As discussed in Section~\ref{sec:hard_soft_at_weak_coupling:soft_interactions:scale}, soft interactions can always be described stochastically, if propagation in the medium is sufficiently long. We quantified this duration as $t \gg \left<\omega^3\right>^2 / \left<\omega^2\right>^3$,
or $S \ll 1$ as defined in Eq.~(\ref{eq:scale}), with $\left<\omega^n\right>$ given by \eqref{eq:rate_moment}.

We can use inelastic interactions
to get an estimate of the length of the medium required to describe soft interactions stochastically. Recall Eq.~(\ref{eq:time_langevin_inel}): it implies that a stochastic description of inelastic interactions with $\omega < T$ requires a time $t \gg 1/[g^4 T]$. For temperatures of a few hundred MeV and a coupling $g\sim 1-2$ encountered in heavy ion collisions,  $1/[g^4 T] < 1$~fm. Under this estimate, it would be reasonable to describe stochastically soft interactions with $\mu \lesssim T$ occurring in a heavy ion collision. Note that this conclusion is based on the estimate of the soft inelastic rate discussed above; should the non-perturbative rate differ significantly from it, it could lead to changes the range of applicability of the Langevin equation. However, we do believe that the current estimates are encouraging.

In what follows, we use $\alpha_s=0.3$ ($g \approx 2$), and first compare a stochastic and a microscopic description of parton energy loss for soft interactions. We use a plasma of length $1$~fm and temperature $T=300$~MeV.

Note that, when the coupling is large, the analytical diffusion coefficients computed perturbatively are not necessarily consistent with numerical values obtained by direct integration of the rates (see Fig.~\ref{fig:qhat_coupling} and surrounding discussion). For what follows, we use the numerical diffusion coefficients in the Langevin part of the hard-soft factorized model.

\subsection{Comparison between diffusion process and collision rate}


\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{plots/soft_collRate_and_langevin_comparison_largecoupling_1fm.pdf}
    \caption{The energy distribution of a $100$~GeV gluon propagating through a $300$~MeV pure glue medium at $\alpha_s = 0.3$. The evolution time is $t = (0.3/\alpha_s )^2=1$~fm. Only soft $1\leftrightarrow 2$ interactions with $\omega < \mu_\omega$ are allowed. Compare with the weak-coupling result from Fig.~\ref{fig:diffusion_collision_rate_small_coupling:inel}.}
    \label{fig:inel_diffusion_collision_rate_large_coupling}
\end{figure*}

\begin{figure*}
    \includegraphics[width=2\columnwidth]{plots/soft_collRate_and_langevin_comparison_elas_large_coupling_numerical_1fm_scale.pdf}
    \onecolumngrid
    \caption{The energy distribution resulting from a $100$~GeV gluon propagating through a $300$~MeV pure glue medium at $\alpha_s = 0.3$. The evolution time is $t = (0.3/\alpha_s )^2=1$~fm. Only soft $2\leftrightarrow 2$ interactions with $\omega < \Lambda$ and $\tilde{q}_\perp < \mu_{\tilde{q}_\perp}$ are allowed. Compare with the weak-coupling result from Fig.~\ref{fig:diffusion_collision_rate_small_coupling:elastic}.}
    \label{fig:elas_diffusion_collision_rate_large_coupling}
\end{figure*}


\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{plots/soft_collRate_and_langevin_comparison_column_elas_large_coupling_200fm.pdf}
    \caption{The energy distribution resulting from a $100$~GeV gluon propagating through a $300$~MeV pure glue medium at $\alpha_s = 0.3$. The evolution time is 200, 50, 20 fm for $\mu_{\tilde{q}_\perp}/T=0.5, 1, 2$; the times were chosen to obtain similarly small values of the skewness parameter $\mathcal{S}$ (\eqref{eq:scale}). Only soft $2\leftrightarrow 2$ interactions with $\omega < \Lambda$ and $\tilde{q}_\perp < \mu_{\tilde{q}_\perp}$ are allowed. }
    \label{fig:diffusion_collision_rate_large_coupling_large_t}
\end{figure*}

As in the weak-coupling case (Section~\ref{sec:weak_coupling_diffusion_vs_rate}), we perform this section's test in the pure glue limit ($N_f=0$). 

We first study the inelastic interactions, and as discussed above, we expect inelastic interactions softer than $\sim T$ to be describable by the Langevin equation in a $1$~fm brick. We show this explicitly in Fig.~\ref{fig:inel_diffusion_collision_rate_large_coupling}.
We show calculations for three different cutoffs $\mu_\omega$,
and we plot the results for the scale $\mathcal{S}$ from \eqref{eq:scale}.\footnote{We verified that the result from \eqref{eq:scale} is close to that of \eqref{eq:scale_inel}. The values we quote are from \eqref{eq:scale}} As expected, agreement between the Langevin approach and the microscopic collision rate approach are best when $\mathcal{S}\ll 1$. In the current setting, agreement is still good for $\mu_\omega = 2 T$, for which $\mathcal{S}=0.33$. This is encouraging evidence that the effect of non-perturbative inelastic interactions ($\mathcal{C}^{1\leftrightarrow 2}_{\textrm{soft}}$) can be treated stochastically in phenomenological applications such as heavy ion collisions.

The equivalent result for soft elastic interactions ($\mathcal{C}^{2\leftrightarrow 2}_{\textrm{soft}}$) is shown in Fig.~\ref{fig:elas_diffusion_collision_rate_large_coupling}. The result is very different. On one hand, the mean energy and width of the parton distribution described with the Langevin equation is almost identical to that described with collision rates.
However their shape are different, especially at smaller values of the cutoffs $\mu_{\tilde{q}_\perp}$. Unlike for the inelastic case and unlike for the elastic case at weak coupling, agreement is improved when the cutoff is larger, which is also reflected in smaller values of the scale $\mathcal{S}$.
The exact reason for this inversion is unknown.
We verified in Fig.~\ref{fig:diffusion_collision_rate_large_coupling_large_t} that longer evolution times do lead to better agreement between the Langevin and the collision rate descriptions, reflected in smaller values of the scale $\mathcal{S}$.
Our conclusion is that soft elastic collision may be more difficult to describe stochastically; it is possible that one needs a larger cutoff $\mu_{\tilde{q}_\perp}$ to describe these elastic interaction stochastically, although more studies will be necessary to confirm this conclusion.
Note, however, that observables which are mainly sensitive to the average energy loss and the width of the parton distribution may tolerate a wider range of soft interactions being described with the Langevin approach. 



\subsection{Parton energy loss at large coupling in a static medium}


\begin{figure*}
    \centering
    \includegraphics[width=2\columnwidth]{plots/cutoff_dependence_large_coupling.pdf}
   
    \caption{The energy distribution of a $100~\textrm{GeV}$ gluon propagating through $300~\textrm{MeV}$ QGP medium ($N_f=3$) at $\alpha_s = 0.3$ with different values of the cutoff. The evolution time is $t=(0.3/\alpha_s)^2=1~\textrm{fm}$. The subplot (a)  only includes $\mathcal{C}^{1\leftrightarrow 2}$ interactions and (b) only includes $\mathcal{C}^{2\leftrightarrow 2}$ interactions. See the weakly-coupled results in Fig.~\ref{fig:cutoff_dependence_small_coupling} for comparison and additional explanations.}
    \label{fig:cutoff_dependence_large_coupling}
\end{figure*}

To close this section, we quantify the cutoff dependence of a $100~\textrm{GeV}$ parton propagating for $1$~fm in a $300$~MeV brick of plasma, with $\alpha_s = 0.3$. This ``brick'' is the same as in the previous section. The soft interactions are described with the Langevin equation, and hard interactions are included as in the full implementation of the hard-soft energy loss model (Section~\ref{sec:factorized_model}). We use $N_f = 3$ in this test. 

We plot the energy distributions with different values of the cutoff in Fig.~\ref{fig:cutoff_dependence_large_coupling}.
In this larger coupling regime, as expected from the results of the previous section, inelastic interactions ($\mathcal{C}^{1\leftrightarrow 2}$) are independent of the cutoff (panel (a)). 
For the elastic case ($\mathcal{C}^{2\leftrightarrow 2}$), the energy distributions with different values of the cutoff are slightly different in the large energy region, although the long tail of the distribution is not affected (panel (b)). 


Note that we also performed a cutoff dependence test on the cutoff $\Lambda$ for $2\leftrightarrow 2$ interactions. We found the energy distribution of a parton propagating in a static medium to be independent of the choice of $\Lambda$, as expected. The result and further discussion can be found in \app{app:lambda_dependence}. 



\section{Application: Energy and fermion-number cascade}

In this section, we use the hard-soft factorized model to study the energy and fermion-number cascade resulting from inelastic interactions between an energetic parton and a thermal medium. This section thus focuses on $\mathcal{C}^{1\leftrightarrow 2}$ (Fig. \ref{fig:reformulation}-b) in the hard-soft factorized model; both the hard splittings and soft splittings are included, with the soft splittings modeled by the Langevin evolution.  The collision kernel $\mathcal{C}^{2\leftrightarrow 2}$ is switched off for this section.


\subsection{Energy cascade of hard gluons}
\label{sec:energy_cascade}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{plots/differential_rate_gg_p1000_Nf0_log.pdf}
    \caption{Comparison between the full leading-order inelastic rate and the deep-LPM regime approximation of the rate from \eqref{eq:branching_rate} with $\hat{q}_{\rm eff}=0.04$~GeV$^3$ for $N_f=0$ ($p=1$~TeV, $T=300$~MeV and $\alpha_s=0.1$).
   
    }
    \label{fig:rate_nf0}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{plots/gg_p1000_different_length_energy_cascade_scaling.pdf}
    \caption{A comparison of the current numerical implementation of QCD kinetics and the analytical approximation of Ref.~\cite{Blaizot:2013hx} for the energy cascade in the pure glue medium for different evolution times. 
    In this test, we only include inelastic $1\leftrightarrow 2$ processes.
   
   
    We use $N_f=0$, $\alpha_s=0.1$, $T=300$~MeV and $p_0=1000$~GeV. }
    \label{fig:energy_scaling}
\end{figure}

When a gluon propagates through a thermal QCD medium, successive medium-induced inelastic radiations result in a gluon cascade. An analytical approximation for the gluon cascade was introduced in Refs.~\cite{Blaizot:2013hx,Blaizot:2015jea}; it was argued that the successive medium-induced quasi-democratic emissions lead to the accumulation of gluons at zero energy and cause a power-law scaling in the small energy region.  We will study this scaling in this section~\cite{[{A pedagogical introduction relating  the AMY kinetic equations to the turbulent cascade presented  here is given in: }]Schlichting:2019abc}.

At leading order, the successive radiations can be assumed to be independent~\cite{Blaizot:2013vha}. In the deep LPM region, where the time scale of the radiation process is much larger than the mean free path between multiple scatterings, the rate per unit time of a gluon with energy $p$ splitting into two gluons with energy fractions $z$ and $1-z$   can be approximated 
as\footnote{Accounting for the identical particles in  the final state, the total rate is $\int_0^{1/2} d\Gamma/dz \,dz$. }~\cite{Blaizot:2013vha, Arnold:2008vd,PhysRevD.79.065025}
\begin{equation}
   \left.\frac{d\Gamma}{dz}\right|_{g\leftrightarrow gg} =\frac{\alpha_s N_c}{\pi}\frac{1}{[z(1-z)]^{3/2}}\sqrt{\frac{\hat{q}_{\rm eff}}{p}}  \, .
\label{eq:branching_rate}
\end{equation}
Here $\hat{q}_{\rm eff}$ is the average transverse momentum broadening of the radiated gluon, and $z=\omega/p$ with $\omega$ the energy of the radiated gluon. 
We have kept only the most singular parts of the splitting 
function at $z \sim 0$.
We  will  treat $\hat q_{\rm eff}$ as a fit parameter, and then relate it to the parameter $\hat q_{\rm soft}^{2\leftrightarrow 2}$ in \eqref{eq:elas_qhat}.


The energy of the initial gluon is $p_0$, and we define $x \equiv \omega/p_0$.  The 
evolution of the gluon spectrum $D(x, \tau) = x(dN/dx)$ is governed by~\cite{Blaizot:2013hx,Schlichting:2019abc}
\begin{multline}
   \frac{\partial D(x, \tau)}{\partial \tau} = \int_0^{1} dz \frac{1}{[z(1-z)]^{3/2}} \\
    \times 
    \left[\sqrt{\frac{z}{x}}D\left(\frac{x}{z}, \tau\right)-\frac{z}{\sqrt{x}}D\left(x, \tau\right)\right], 
    \label{eq:FP_evolution}
 \end{multline}
where
\begin{equation}
   \tau \equiv \frac{\alpha_sN_c}{\pi}\sqrt{\frac{\hat{q}_{\rm eff}}{p_0}}t \, ,
\end{equation}
and $t$ is the evolution time of the gluon. 

The exact solution for \eqref{eq:FP_evolution} can be calculated via Laplace transform:
\begin{equation}
    D_0(x, \tau)=\frac{\tau}{\sqrt{x}(1-x)^{3/2}}e^{-\pi[\tau^2/(1-x)]}. 
    \label{eq:spectrum}
\end{equation}
As remarked in Ref.~\cite{Blaizot:2013hx}, this power-law gluon spectrum Eq. (\ref{eq:spectrum}) scales as $1/\sqrt{x}$ in the small-x region.  

In order to compare  with \eqref{eq:spectrum}, we first determine the approximate value of $\hat{q}_{\rm eff}$ to use in the simplified rate \eqref{eq:branching_rate}; this value also enters Eqs.~(\ref{eq:FP_evolution}-\ref{eq:spectrum}).
We fix $\hat{q}_{\rm eff}$ by comparing \eqref{eq:branching_rate} with the full leading-order inelastic rate, as shown in Fig.~\ref{fig:rate_nf0}. With parameters given in Fig.~\ref{fig:rate_nf0}, we find  $\hat{q}_{\rm eff}\simeq 0.04$~GeV$^3$ at $\omega/p\simeq 10^{-2}$. We will use this value of $\hat q_{\rm eff}$ in our analysis of the cascade below.

It should be emphasized that \eqref{eq:branching_rate} is an approximation to the full inelastic rates corresponding to \eqref{eq:1to2}. Indeed, a leading-log analysis of the full rates at small $z$ in the deep LPM regime shows that~\cite{Arnold:2008vd}:
\begin{equation}
   \hat q_{\rm eff} = \hat q_{\rm soft}^{2\leftrightarrow 2} (\mu_\perp^2) \, ,
   \label{eq:qeff_est}
\end{equation}
where  $\hat{q}_{\rm soft}^{2\leftrightarrow 2}$ is given in \eqref{eq:elas_qhat}, and  $\mu_\perp^2=C_0\sqrt{2 \omega \hat q_{\rm eff} }$ with $C_0 \sim 1$.
The cutoff $\mu_\perp^2$ scales with 
the accumulated
transverse momentum of the radiated gluon over its formation time. A next-to-leading logarithmic analysis fixes the coefficient $C_0$~\cite{Arnold:2008zu}: 
\begin{equation}
\mu_\perp^2 = C_0 \sqrt{2 \omega \hat{q}_{\rm eff}}\, ,   \qquad C_0 = 2e^{2 -\gamma_E + \pi/4} \, .
\label{eq:muperp_est}
\end{equation}

For $N_f=0$, $p=1$~TeV, $T=300$~MeV, $\alpha_s=0.1$ (same as in \Fig{fig:rate_nf0}), and using $\omega/p=10^{-2}$, we can solve Eqs.~(\ref{eq:qeff_est}-\ref{eq:muperp_est}) numerically. We find $\hat{q}_{\rm eff}\approx 0.052$, which as expected is close to the value we found in Fig.~\ref{fig:rate_nf0}.



We next perform the gluon cascade in a pure-glue medium ($N_f=0$) using the hard-soft factorized model, i.e we  include both soft inelastic interactions described with the Langevin equation, and rate-based hard inelastic interactions which dominate this test. 
In Fig. \ref{fig:energy_scaling}, we compare this numerical result calculated by the current model with the analytical spectrum in Eq. (\ref{eq:spectrum}). We find that the numerical solution for the medium-induced cascade is reasonably well-described by the approximate analytic solution. In particular, the power 
law behavior, $dN/dx \propto x^{-3/2}$, is nicely captured by this solution.


\subsection{Fermion-number cascade of gluons and quarks}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{plots/fermion_number_ratio_10TeV_QGP.pdf}
    \caption{The fermion number cascade of the numerical implementation in the QGP medium with different evolution times. In this test, we only include inelastic interactions ($\mathcal{C}^{1\leftrightarrow 2}$). We use $N_f=3$, $\alpha_s=0.3$, $T=300$~MeV and $p_0=10$~TeV. The black horizontal line reflects the expected limiting value of $D_S/(2 N_f D_g) \approx 0.07$ (\eqref{eq:fermion_ratio}).}
    \label{fig:fermion_cascade}
\end{figure}

The fermion-number cascade was investigated in Ref.~\cite{Mehtar-Tani:2018zba}. Given the power-law scaling in the small energy region, at small $x$, we can write the power-law spectrum of quarks and gluons as
\begin{equation}
    \begin{split}
        D_g &\equiv x \frac{dN_g}{dx} = \frac{G}{\sqrt{x}},  \\
        D_s &\equiv \sum_{i = 1}^{N_F}\left(D_{q_i}+D_{\bar{q}_i}\right) = \frac{Q}{\sqrt{x}}. 
    \end{split}
\end{equation}

As derived in Ref.~\cite{Mehtar-Tani:2018zba}, the quark-to-gluon ratio of the soft radiated partons is determined by the transformation rate between gluons and fermions. We have
\begin{equation}
    \frac{Q}{2N_fG} = \frac{1}{2N_f}\frac{\int_0^1 dz z \mathcal{K}_{qg}(z)}{\int_0^1 dz z \mathcal{K}_{gq}(z)} \approx 0.07, 
    \label{eq:fermion_ratio}
\end{equation}
where $\mathcal{K}_{qg}$ is the splitting function of $g \rightarrow q\bar{q}$, and $\mathcal{K}_{gq}$ is the splitting function of $q \rightarrow gq$. 

To test the quark-to-gluon ratio in the hard-soft factorized model, we numerically simulate the evolution of a gluon or a quark propagating through a static QGP medium ($N_f=3$) using the full leading order inelastic rate. 
We perform the calculation for both an energetic gluon and an energetic quark with an initial energy of $10$~TeV.
The result is shown in Fig. \ref{fig:fermion_cascade}; we find that it converges to the universal quark-to-gluon ratio when using the full QCD rates. 


\section{Summary and outlook}

This work introduces a new formulation of parton energy loss where soft and hard interactions with the underlying plasma are factorized and treated separately.
The factorization is performed with cutoffs based on the momentum transfer of the interactions.
Rare hard interactions are considered as independent successive interactions, and solved with collision rates (Sections~\ref{sec:factorized_model:inelastic_hard} and \ref{sec:factorized_model:elastic_hard});
the larger momentum exchange with the medium make them more likely to be amenable to a perturbative description.
On the other hand, frequent soft interactions are treated stochastically using a Langevin evolution with drag and diffusion coefficients encoding the effect of these soft interactions (Section~\ref{sec:hard_soft_theory:reformulation:soft}); non-perturbative effects can thus be absorbed in these transport coefficients. 

Our numerical implementation of this model (Section~\ref{sec:tequila_weak_coupling}) shows that this factorization works well in the weakly-coupled regime where the theory was derived~\cite{Ghiglieri:2015ala}.
In fact, by revisiting the conditions under which the Langevin equation can describe the Boltzmann equation (Section~\ref{sec:hard_soft_at_weak_coupling:soft_interactions:scale}), we extended the region of phase space (``cutoffs'') that can be described stochastically. We used the dimensionless scale $\mathcal{S}$ (\eqref{eq:scale}) to quantify the length of a plasma necessary for soft collisions to be describable with the Langevin equation. Our numerical tests showed that this scale works very well in practice.

Because the scale $\mathcal{S}$ is a property of the Boltzmann equation and not a perturbative concept, we used it to extend our discussion of parton energy loss beyond the perturbative regime. We estimated that inelastic collisions resulting in parton energy loss of order $T$ could be described stochastically in a QCD plasma of size $\sim 1$~fm (Section~\ref{sec:large_coupling}). Given that inelastic interactions dominate parton energy loss for high-energy partons, this supports the applicability of the present energy loss model in heavy ion collisions.


This work paves the way to systematic phenomenological constraints on the soft transport coefficients of light partons.
The key strength of our approach is that perturbative parton energy loss calculations are still being used for harder interactions --- the regions of phase space where they are most likely to hold.
Conversely, the interactions most sensitive to non-perturbative effects --- soft interactions --- are encoded in simple transport coefficients which can be constrained by comparison with measurements.
A stochastic description of soft collisions can also be very efficient numerically, as a large number of soft interactions can be absorbed in the transport coefficients.
These phenomenologically-constrained transport coefficients can eventually be compared with  lattice results~(e.g. Ref.~\cite{Moore:2019lgw}).
A similar program is already being pursued for the energy loss of heavy quarks~\cite{Ke:2018tsh}; studies of light parton energy loss with a model that includes many features of soft-hard factorization, are also ongoing~\cite{Ke:2020clc}.

Future generalization of this framework includes improving the treatment of the radiation angle of collinear radiation, and the inclusion of a running coupling constant and of next-to-leading order effects; these additions will increase the type of observables that can be studied with this model. The inclusion of finite-size effects in this formalism will also be an important addition. These additions will be able to build on Ref.~\cite{Ke:2020clc} and other works.


\begin{acknowledgments}
We thank Weiyao Ke for his invaluable help in the early stage of this project, and Jacopo Ghiglieri for generously sharing notes on the splitting approximation that formed the basis of this work's discussion on the topic.
We thank Sangyong Jeon, Chanwook Park, Abhijit Majumder and the other members of the JETSCAPE Collaboration for discussions regarding MARTINI and MATTER, and for their support with the JETSCAPE framework.
We thank Yingru Xu for valuable discussions regarding the Langevin equation and its numerical implementation.
This work was supported by the U.S. Department of Energy Grant no. DE-FG02-05ER41367 (SAB, JFP and TD) and DE-FG-02-08ER41450 (DT).
TD is also supported by NSF grant OAC-1550225. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.
\end{acknowledgments}

