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


\section{Introduction and Motivation}
The concept of aerosol mixing state describes how different aerosol
chemical species are distributed among the particles in a population
\citep{RIEMERETAL19}. A completely ``externally mixed'' population
contains only one species per particle, while a completely
``internally mixed'' population contains particles with the same
composition as the bulk population. Many different intermediate mixing
states are possible between those two extremes, and these commonly
occur in the ambient atmosphere. Aerosol mixing state influences the
particles' properties such as hygroscopicity
\citep{FIERCEETAL17,HOLMGRENETAL14}, optical properties
\citep{FIERCEETAL16,JACOBSON01}, cloud condensation nuclei activity
\citep{CHINGETAL12,WANGETAL10}, ice nucleation potential
\citep{KNOPFETAL18}, the aerosols' lifetime \citep{KOCHETAL09} in the
atmosphere, and their deposition in the human respiratory system
\citep{CHINGKAJINO18}. However, Earth System models (ESMs) or regional
climate models usually hold simplified mixing state assumptions. These
can influence how accurately physicochemical properties of aerosols
are predicted \citep{CHINGKAJINO18, FIERCEETAL17}, and thereby limit
the accuracy of estimating aerosol impacts on climate and health.

\citet{RIEMERWEST13} proposed an entropy-based diversity metric to
quantify aerosol mixing state. This metric, the mixing state index
($\chi$), has a range from 0\% (a fully externally mixed population)
to 100\% (a fully internally mixed population). The index $\chi$
has been applied to both modeling \citep{ZHUETAL16} and experimental
\citep{HEALYETAL14,OBRIENETAL15,FRAUNDETAL17,BONDYETAL18,YEETAL18}
work. For example, $\chi$ has been used for
quantifying errors in CCN concentration \citep{CHINGETAL17} prediction,
as well as the prediction of soot particles depositing in the human
respiratory system \citep{CHINGKAJINO18}.

Knowing the global distribution of aerosol mixing state index is
desirable, as it can be used for error quantification at the large
scale. But deriving this distribution directly by using a benchmarking
particle-resolved model \citep{RIEMERETAL09} at the global scale would
be computationally prohibitive over the time scales of interest. To
overcome this limitation, supervised-learning emulators were developed
for predicting the spatial distribution aerosol mixing state
across the globe \citep{HUGHESETAL18,ZHENGETAL20}. However, the
spatial delimitation based on aerosol mixing state indices remains
unclear, which hinders further spatiotemporal analysis.

The research questions for this paper are: (1) Is it possible to
regionalize the global mixing state indices, i.e., define aerosol
mixing state zones, similar to distinct climate zones
\citep{FOVELLFOVELL93,MITCHELL76}? (2) What are spatiotemporal
patterns of the regionalized mixing state indices? This paper
describes an effort to regionalize the global mixing state indices
using a simple but effective unsupervised learning approach.

\section{Methods}
\subsection{Mixing state indices definition}

The mixing state index $\chi$ \citep{RIEMERWEST13} is given by the affine
ratio of the average particle species diversity ($D_{\alpha}$) and bulk
population species diversity ($D_{\gamma}$) as
${\chi} = \frac{D_{\alpha}-1}{D_{\gamma}-1}$.
The diversities $D_{\alpha}$ and $D_{\gamma}$ are calculated based on
per-particle species mass fractions and mass fraction of particles as
described in detail in \citet{RIEMERWEST13}.

This study focused on submicron aerosols because they are the dominant
category of particles for radiation interactions and the provision of
cloud condensation nuclei. We defined the mixing state indices in
three different ways,
namely based on chemical species abundance ($\chi_\textrm{a}$), based
on the mixing of absorbing and non-absorbing species
($\chi_\textrm{o}$), and based on the mixing of hygroscopic and
non-hygroscopic species ($\chi_\textrm{h}$). The index
$\chi_\textrm{a}$ was defined based on all the model aerosol
species. A lower value for $\chi_\textrm{a}$ indicates that the
species tend to be present in separate particles. For
$\chi_\textrm{o}$, we considered two surrogate species, black carbon
(absorbing) and all other aerosol species grouped together (assumed to
be non-absorbing). Thus, a lower value for $\chi_\textrm{o}$ indicates
the absorbing species black carbon and the sum of all other
(non-absorbing) species are more externally mixed. Similarly,
$\chi_\textrm{h}$ was also calculated from two surrogate species. We
combined black carbon, dust, and primary organic matter as one
surrogate species, given their comparatively lower hygroscopicities,
and salt, secondary organics, and sulfate as the other surrogate
species. Correspondingly, a lower value for $\chi_\textrm{h}$
represents the case where hygroscopic and non-hygroscopic species
tended to be present in separate particles. Note that mixing state
indices defined for different purposes capture different aspects of
the overall mixing state, and they are uncorrelated
\citep{ZHENGETAL20}.

\subsection{Data}
We used the supervised-learning surrogate models developed by
\citet{ZHENGETAL20} to predict the mixing state indices from ESM
(Earth System Model) output. The surrogate models were trained on an
ensemble of particle-resolved model simulations generated by the
particle-resolved model PartMC-MOSAIC
\citep{RIEMERETAL09,ZAVERIETAL08}, with the training labels for $\chi$
calculated directly from the particle-resolved data. The strategy to
generate the training and testing data was to vary 45 input parameters
for the PartMC-MOSAIC model scenarios, including primary emissions of
different aerosol types, primary emissions of gas phase species, and
meteorological parameters. A Latin Hypercube sampling approach was
employed to sample the parameter space efficiently. The surrogate
models were trained by using XGBoost (eXtreme Gradient Boosting,
\citet{CHENGUESTRIN16}), and can be expressed as:
\begin{equation}
\chi_S = f_S(A,G,E).
\end{equation}
Where $\chi_S$ is the mixing state index ($\chi_\textrm{a}$,
$\chi_\textrm{o}$, or $\chi_\textrm{h}$) and $f_S$ denotes
the surrogate function for the corresponding mixing state index.
The variants $A$ (aerosol), $G$ (gas), and $E$
(environment) represent the features (variables) from the PartMC
simulations that are used for predicting the mixing state index. These variables are
also available from the ESM.  A detailed definition of the features for the
surrogate models is given in \citet{ZHENGETAL20}.

We used the simulations from the Community Earth System Model Version
2 \citep[CESM2 version 2.1.0;][]{DANABASOGLUETAL20} with MAM4
\citep{LIUETAL16} to provide variables at the global
scale. Specifically, we used the component set ``FHIST'' for the
global simulation configuration. This component set represents a
typical historical simulation using an active atmosphere and land with
prescribed sea surface temperatures and sea-ice extent, and a 1~degree
finite volume dycore with the forcing data available from 1979 to
2015. We ran the model for the year 2011 with 6~years (2005-2010)
spinup. The simulation was conducted at a resolution of 0.9$^\circ$
latitude by 1.25$^\circ$ longitude along with emission inventories
from CMIP6 emissions \citep{EMMONSETAL20}. We stored the instantaneous
outputs every three hours in the simulation, which yielded
2920~timestamps for each surface-layer grid cell for the entire year
of simulation time.

With the surrogate models and the ESM simulation, we predicted the
mixing state indices for each grid cell at each timestamp (i.e. every
3~hours). For each grid cell, the mixing state indices were averaged by
month.  Therefore, each grid cell has a vector of stacked mixing state
indices with a length of 36 (12 months $\times$ 3 mixing state types),
represented as $\boldsymbol{x} =
(\chi_{\textrm{a},1},\dots,\chi_{\textrm{a},12},\chi_{\textrm{o},1},\dots,\chi_{\textrm{o},12},\chi_{\textrm{h},1},\dots,\chi_{\textrm{h},12})$.

\subsection{Regionalization strategy}
Here we interchangeably use the terms ``region'' and ``cluster'' to
refer to a group of grid cells. We used the $k$-means unsupervised
clustering algorithm to partition $n$ grid cells into $k$ clusters,
minimizing the within-cluster variances. The $k$-means method has
been applied to environmental sciences for ecoregion delineation
\citep{KUMARETAL11}, environmental risk zoning of the chemical
industrial area \citep{SHIZENG14}, and clustering haze trajectory of
peatland fires \citep{KHAIRATETAL17}, among other applications.

In our case, given $n=55\,296$ (192 latitude $\times$ 288 longitude) grid cells
$(\boldsymbol{x}_1, \boldsymbol{x}_2, \dots, \boldsymbol{x}_{55296})$, where
each grid cell contains a 36-dimensional vector (three mixing state indices over
12 months), $k$-means clustering aims to partition these grid cells into $k$
($k \le n$) sets $\boldsymbol{S} = \{S_1,S_2,\dots,S_k\}$ so as to
minimize the within-cluster variance. The objective is to find:
\begin{equation}
\mathop{\arg\min}_{\boldsymbol{S}} \ \ \ \sum_{i=1}^{k}\sum_{\boldsymbol{x}\in S_i}=\Vert{\boldsymbol{x}}-\boldsymbol{\mu}_i\Vert^2
\end{equation}
where $\boldsymbol{\mu}_i$ is the mean of the points in $S_i$.

The $k$-means clustering method requires the a priori specification of
the number of clusters ($k$). As the value of $k$ increases, the
variance $V_k$ asymptotically approaches 0 until $k$ equals $n$. Here
we used the relative difference in variance between $k-1$ and $k$ to
determine the number of clusters. Starting from $k=2$, we computed the
criteria of relative difference $V_{k-1} - V_{k} \le 0.05V_{k-1}$ to
determine the cluster number $k$. In this study, since
$V_{10} - V_{11} \le 0.05V_{10}$, we selected $k=11$
clusters.

\section{Results}

The regional clusters of mixing state are shown in
Figure~\ref{fig:map_clusters}. Unlike with supervised learning, the
performance of unsupervised learning are challenging to evaluate since
there is no ground truth.  However, the clustering process of this
study was merely based the information of mixing state indices,
without the direct guidance such as aerosol species distribution or
other explicit spatial information. Therefore a success indicator of
the spatiotemporal clustering is whether spatially-contiguous regions
emerge from this procedure. The clusters in
Figure~\ref{fig:map_clusters} exhibit notable regional patterns among
the eleven clusters (mixing state zones), suggesting that this
approach has indeed identified meaningful clusters. Each cluster has a
unique spatiotemporal distribution of three mixing state indices, and
captures the spatiotemporal variations that could not be detected by
the global overall averages. For example, ocean regions (e.g.,
clusters 8, 9, 10 in the Southern Hemisphere) tend to display a more
pronounced seasonal cycle than land regions (e.g., cluster 7) and
overall averages.

Clusters~1 (the Arctic) and~11 (the Antarctic) were characterized by
high mixing state indices (with monthly means of at least 67\%)
throughout the year. The high values of all three mixing state indices
means that the different aerosol species are rather internally mixed
in these regions. Clusters~5 (tropical oceans), 8 (the South China
Sea, Indian Ocean and the Pacific Ocean in the northern hemisphere
near the equator, and the northern part of the Southern Ocean) and~9
(the southern part of the Southern Ocean) were ocean areas. Similar
bulk aerosol composition was predicted in these regions, however the
annual cycle of mixing state varies more for clusters~8 and~9 than
for~5. Clusters~2, 3, and~4 covered oceans at mid and high latitudes
in the Northern Hemisphere, and lands including Northern and Southern
America, Europe, Middle East, Southeast Asia, as well as
Australia. Clusters~6 (the region
across the Atlantic from Western and Northern Africa to the northeast
coast of South Africa, as well as part of the Indian Ocean) and~7
(Middle, Eastern, and Southern Africa, as well as Asia) were
characterized by their higher values in $\chi_\textrm{o}$ compared to
$\chi_\textrm{a}$ and $\chi_\textrm{h}$, meaning that black carbon and
non-absorbing species are comparatively internally mixed.

\begin{figure}
  \centering
  \includegraphics[width=\textwidth,trim=50 0 0 0]{./mixing_state_zones.pdf}
  \caption{Aerosol mixing state zones based on unsupervised
    regionalization. Also shown are monthly average mixing state
    indices and annual averages of aerosol species mass fractions
    (bulk composition).}\label{fig:map_clusters}
\end{figure}

\section{Conclusions}
In this paper, we developed a simple but effective unsupervised
learning approach to regionalize global aerosol mixing state
indices. We used the monthly averages of the spatially-varying aerosol
mixing state indices as input data. Each grid cell was then assigned
to a region using $k$-means clustering without explicit spatial
information as input. This approach resulted in eleven distinct
regions over the world with specific spatial aggregation
patterns. Each region exhibited a unique distribution of mixing state
indices, suggesting that the unsupervised regionalization approach had
identified meaningful regions. To the best of our knowledge, this is
the first study to define aerosol mixing state zones, which we
anticipate will be helpful for future studies of the global aerosol
burden.

\begin{ack}
We would like to acknowledge high-performance computing support from
Cheyenne (\url{doi:10.5065/D6RX99HX}) provided by NCAR's Computational and
Information Systems Laboratory, sponsored by the National Science
Foundation. The CESM project is supported primarily by the National
Science Foundation.
This research used resources of the Oak Ridge Leadership Computing
Facility, which is a DOE Office of Science User Facility supported
under Contract DE-AC05-00OR22725. This research was supported in part
by an appointment to the Oak Ridge National Laboratory ASTRO Program,
sponsored by the U.S. Department of Energy and administered by the Oak
Ridge Institute for Science and Education. We also acknowledge funding
from DOE grant DE-SC0019192 and NSF grant AGS-1254428.
This research is part of the Blue Waters sustained-petascale computing
project, which is supported by the National Science Foundation (awards
OCI-0725070 and ACI-1238993) the State of Illinois, and as of
December, 2019, the National Geospatial-Intelligence Agency. Blue
Waters is a joint effort of the University of Illinois at
Urbana-Champaign and its National Center for Supercomputing
Applications. Louisa Emmons is thanked for thoughful comments on the
CESM2 simulations and the manuscript. We thank AWS for providing AWS
Cloud Credits for Research.
\end{ack}

\clearpage


\section{Summary}

WI Fast Plants (\url{https://fastplants.org/}) \cite{williams1986rapid} were developed as a research tool at the University of Wisconsin-Madison and have been used by K-16 teachers around the world for nearly 30 years as an educational model-organism. As the name suggest, WI Fast Plants have a short time from planting to flowering (about 2 weeks), and thus, these plants are used in classrooms to engage learners of all ages and grade levels into the investigation of plant life cycles, the growth of flowers, the role of environmental factors on plants, the energy transformation of plants, the genetics of hybrids among many other topics.

Thousands of students from elementary school to college level grow WI Fast Plants in class or at home, and collect data from their experiments, yet they do not have any user-friendly tool to visualize or analyze these data. Existing Data Science platforms such as CODAP \cite{codap} are either too complicated to use or expensive. Furthermore, existing interfaces are not tailored to WI Fast Plants data or educational objectives, and lack the flexibility to evolve alongside the lesson options provided by WI Fast Plants.

\section{Statement of need}

WI Fast Stats (\url{https://wi-fast-stats.wid.wisc.edu/}) is an integrated animated web page which serves as a medium to a collection of R-developed web apps that provide Data Visualization and Data Analysis tools for WI Fast Plants data.
Each web app corresponds to a educational unit linked to a specific WI Fast Plants webinar and it serves two main functions: 1) K-16 teachers attending the WI Fast Plants webinar will learn how to design Data Science exercises through the web app for their students and 2) students can use the web app independently to learn about visualization and analysis via the publicly available sample datasets and educational materials.

In this sense, WI Fast Stats is designed to be used by both K-16 teachers and by
students from elementary school to college. A preliminary version of a web app has already been used in a 2020 WI Fast Plants webinar with over 300 high-school teachers across the U.S. in attendance. Based on the comments at the webinar and follow-up messages, the attendees were deeply impressed by the potential of the web app to revolutionize the manner in which they currently teach Data Science in the classroom.

WI Fast Stats is a user-friendly easy-to-use interface that will render Data Science accessible to teachers and students without strong programming or mathematical background. Because WI Fast Stats was also created and is maintained at the University of Wisconsin-Madison, it is the perfect companion to the already widely successful educational tool of WI Fast Plants. WI Fast Stats will evolve alongside WI Fast Plants by fulfilling the Data Science needs of the thousands of teachers and students around the world who currently use WI Fast Plants to learn the complexities of plants and the biological world.

\section{Description of the website}

WI Fast Stats (\url{https://wi-fast-stats.wid.wisc.edu/}) is an integrated animated web page which hosts a collection of web applications and relevant webinars for WI Fast Plants. The web page contains six modules including Home, About, Webapps, Webinars, Source Code, and FAQ (Figures \ref{figweb1} and \ref{figweb2}). The Home page utilizes \texttt{animate.css} to construct the animation with three scrolling images retrieved from WI Fast Plants. \texttt{Font-awesome.css} offers a fantastic mode to show buttons representing the links of our web applications as well as the Github source code. We also built a timeline to store all the webinar-related information.

The web page and accompanying web apps are all open-source with the code stored in the GitHub repository: \url{https://github.com/crsl4/fast-stats}.

\begin{figure}
    \centering
    \includegraphics[scale=0.25]{website1.png}
    \caption{WI Fast Stats website home page: \url{https://wi-fast-stats.wid.wisc.edu/}.}
    \label{figweb1}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[scale=0.25]{website2.png}
    \caption{WI Fast Stats website comprises different web apps each corresponding to a given WI Fast Plants webinar.}
    \label{figweb2}
\end{figure}

\section{Description of the web apps}

The web apps within WI Fast Stats provide an interactive and easy-to-use platform for data visualization and data analysis for data collected in accordance to the WI Fast Plants webinars and educational materials.
The web apps are based on the \texttt{R shiny} package and contain three main modules: Data Summary, Data Visualization, and Data Analysis which are built with the \texttt{ShinyDashboard} framework (Figures \ref{figwebapp1} and \ref{figwebapp2}).

The Data Upload section allows the users to upload their own collected data based on their fast plants experiments, or to utilize the already loaded sample dataset which illustrates the same educational outcomes intended in the webinar without having to run the experiments.
In addition, this section also provides a summary button to show only first rows of the dataset or the whole dataset.

The Data Visualization section allows the users to create five plots: Mosaic plot, Scatter plot, Box plot, Violin plot, and Density plot. It utilizes the \texttt{ggplot2} library to declaratively create graphs based on particular group variables and quantity variables. Furthermore, the \texttt{Plotly} library extends the features of the plot by adding \texttt{Lasso Select}, \texttt{autoscale}, and \texttt{data-toggle} tooltip. Finally, color palettes, transparency, point size, and point shape are available options to improve the overall appearance of the graph.

The Data Analysis section (not present in web apps tailored to middle
school students) allows the user to perform statistical tests like the
chi-square test and the t test to compare the characteristics of the plants under different environmental or experimental settings.

Finally, the web apps maintain a validator system detecting any illegal actions done by users and providing meaningful error messages. The website and web apps are accompanied by a specialized Google user group (\texttt{wi-fast-stats}) for general questions.

\begin{figure}
    \centering
    \includegraphics[scale=0.25]{webapp1.png}
    \caption{WI Fast Stats web app corresponding to the WI Fast Plants webinar on the selection of polycot plants: \url{https://wi-fast-stats.wid.wisc.edu/cotyledon/}.}
    \label{figwebapp1}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[scale=0.25]{webapp2.png}
    \caption{WI Fast Stats web app corresponding to the WI Fast Plants webinar on the effect of the ecosystem on the plants.: \url{https://wi-fast-stats.wid.wisc.edu/ecosystem/}.}
    \label{figwebapp2}
\end{figure}

\section{Future work}

WI Fast Stats will be continuously evolving to provide data visualization and data analysis capabilities for the ever-growing needs of the WI Fast Plants community.

\section*{Acknowledgements}
This work was supported by the Department of Energy [DE-SC0021016 to CSL]. We thank Hedi Baxter Lauffer and everybody at WI Fast Plants for inviting us to work with them on the creation of these Data Science educational open-source tools. Finally, we acknowledge the work in \cite{Hotaling2020} which helped us improve the scientific writing of this manuscript.

\bibliographystyle{unsrt}