# Parton shower algorithm with saturation effect

Yu Shi  Key Laboratory of Particle Physics and Particle Irradiation (MOE),
Institute of Frontier and Interdisciplinary Science, Shandong University,
Qingdao, Shandong 266237, China  Shu-Yi Wei  Key Laboratory of Particle
Physics and Particle Irradiation (MOE), Institute of Frontier and
Interdisciplinary Science, Shandong University, Qingdao, Shandong 266237,
China  Jian Zhou  Key Laboratory of Particle Physics and Particle Irradiation
(MOE), Institute of Frontier and Interdisciplinary Science, Shandong
University, Qingdao, Shandong 266237, China

###### Abstract

We extend the previously developed small $x$ parton shower algorithm to
include the kinematic constraint effect and $k_{t}$ resummation effect. This
work enables the Monte Carlo generator to simultaneously resum large $k_{t}$
and small $x$ logarithms in the saturation regime for the first time. It is an
important step towards simulating processes involving multiple well separated
hard scales, such as di-jet production in eA collisions at EIC.

##  I  Introduction

The study of dense gluonic matter at small $x$ inside a large nucleus and
nucleon has been and continues to be an important frontier of high-energy
nuclear physics. It is also one of the main objectives of the physics program
of the future Electron-Ion Collider (EIC) Accardi  _et al._ (2016); Abdul
Khalek  _et al._ (2022). Tremendous theoretical efforts have been made to
search for smoking gun evidence of saturation. To this end, hard scattering
processes in eA collisions at EIC are expected to deliver crucial messages
about how saturation emerges from strongly interacting gluonic matter. A Monte
Carlo event generator that incorporates saturation effects could play an
essential role in fully harnessing the potential of future experimental data
taken from EIC.

As the core of general purpose Monte Carlo event generators, parton showers
describe successive radiations from highly-energetic partons that participate
in the hard scattering process. While most parton branching algorithms Ellis
_et al._ (2011); Buckley  _et al._ (2011); Höche (2015); Byer  _et al._ (2022)
are based on the soft and collinear approximation which effectively resums the
Dokshitzer-Gribov-Levin-Altarelli- Parisi (DGLAP) Altarelli and Parisi (1977)
like logarithm to all orders, only a few parton shower generators Marchesini
and Webber (1991, 1992); Jung and Salam (2001); Jung  _et al._ (2010a);
Andersen and Smillie (2010a, b); Jung  _et al._ (2010b); Andersen and Smillie
(2011); Lipatov  _et al._ (2020); van Hameren (2018); Hautmann  _et al._
(2022) have been developed to describe small $x$ processes by simulating semi-
hard emissions which give rise to the logarithm of the type $\ln(1/x)$ Kuraev
_et al._ (1977); Balitsky and Lipatov (1978). Among these generators, the
Cascade Jung and Salam (2001); Jung  _et al._ (2010a) that is built on the
Catani-Ciafaloni-Fiorani-Marchesini (CCFM) evolution equation Ciafaloni
(1988); Catani  _et al._ (1990a, b); Marchesini (1995) is the most widely used
in the phenomenology studies (see for recent examples Lipatov  _et al._
(2023); Lipatov and Malyshev (2023)). However, none of the aforementioned
parton showers takes into account the gluon recombination process occurs in
the dense target.

The first attempt to include saturation effect in the parton shower is
presented in Ref. Shi  _et al._ (2023) where both the forward and the backward
evolution schemes have been presented. The underlying parton branching
equation employed in our formulation is the folded Gribov-Levin-Ryskin (GLR)
equation Gribov  _et al._ (1983). Although the GLR equation is somewhat
outdated compared to modern treatments of small $x$ evolution Balitsky (1996);
Kovchegov (1999); Jalilian-Marian  _et al._ (1997, 1998); Iancu  _et al._
(2001); Ferreiro  _et al._ (2002), it is sufficient for simulating events in
eA collisions at EIC energy. This is because the gluon density probed at EIC
is not high enough for the triple pomeron vertex to dominate the gluon fusion
process. In the previous work Shi  _et al._ (2023), we performed a consistent
check by comparing the transverse momentum distribution of exchanged gluons
reconstructed from the parton shower generator with numerical solutions of the
GLR equation. A full agreement between these two results was reached. The
running coupling effect was also implemented in our Monte Carlo simulation.

In the present work, we improve this parton branching algorithm by imposing
the kinematic constraint arising from the requirement that the offshellness of
$t$ channel gluon should be dominated by its transverse momentum squared
Kwiecinski  _et al._ (1996, 1997); Deak  _et al._ (2019). Though it is
formally a sub-leading logarithm contribution, the kinematic constraint effect
is known to significantly slow down the evolution speed. It is thus a
necessary component of the Monte Carlo generator for any practical
phenomenological studies. Actually, the angular ordering of soft emissions is
automatically imposed once the kinematic constraint is applied since the
angular ordering constraint is weaker than the latter Kwiecinski  _et al._
(1996) in the small $x$ limit. The coherent branching effect is thus
effectively included in the parton shower. On the other hand, for the case of
hard scattering processes involving multiple well-separated hard scales, like
di-jet production in eA collisions, the transverse momentumn dependent (TMD)
type large logarithm $\alpha_{s}\ln^{2}\left(Q^{2}/k_{\perp}^{2}\right)$ and
small $x$ logarithm $\alpha_{s}\ln\left(1/x\right)$ need to be simultaneously
resummed. Such a joint resummation formalism has been established in a series
of publications Zhou (2016); Xiao  _et al._ (2017); Zhou (2019). Another main
objective of this work is to implement the joint resummation in the Monte
Carlo simulation.

The rest of the paper is organized as follows. In Sec. II, we discuss how to
integrate the kinematic constraint effect into the parton shower algorithm.
The formulations of both forward and backward evolution are presented. In Sec.
III, the implementation of the joint resummation in the algorithm is
discussed. Our starting point is the Sudakov factor derived from a folded
version of the Collins-Soper (CS) and the renormalization group equation. It
is shown that the $k_{\perp}$ distribution reconstructed from the parton
shower is identical to the numerical and analytical results obtained from the
CS equation and renormalization group equation. The paper is summarized in
Sec. IV.

##  II The kinematic constraint

In our previous work Shi  _et al._ (2023), we have developed a Monte Carlo
method to simulate the parton shower at small $x$ based on the GLR evolution
equation Gribov  _et al._ (1983). Our formulation only takes into account the
summation of the leading logarithm $\ln\left(1/x\right)$ contribution which is
known to result in too rapid growth of gluon number density towards small $x$
region. From a phenomenological point of view, it is crucial to go beyond the
leading logarithm accuracy and include the various sub-leading logarithm
contributions Kwiecinski  _et al._ (1996, 1997); Deak  _et al._ (2019); Liu
_et al._ (2022); Avsar  _et al._ (2011); Iancu  _et al._ (2015a, b); Lappi and
Mäntysaari (2016); Ducloué  _et al._ (2019, 2020); Zheng and Zhou (2019),
among which the kinematic constraint effect Kwiecinski  _et al._ (1996, 1997);
Deak  _et al._ (2019); Liu  _et al._ (2022) is a particularly interesting one.
The kinematic constraint is required for the validity of the BFKL/GLR equation
at small $x$. The constraint is needed to ensure that the virtuality of the
gluons along the chain is controlled by the transverse momenta. The
implementation of the kinematic constraint can significantly slow down the
small $x$ evolution and thus lead to a better description of relevant
phenomenology. Note that the angular ordering of the gluon emissions is
automatically satisfied once the kinematic constraint is imposed in the small
$x$ limit. The coherent branching effect is thus effectively achieved
following the steps outlined below.

The starting point of the Monte Carlo implementation for such an effect is the
folded GLR equation with the kinematic constraint. Following the arguments
made in Refs. Kwiecinski  _et al._ (1996); Deak  _et al._ (2019), the
transverse momentum square of the radiated gluon $l_{\perp}^{2}$ must be
smaller than $\frac{1-z}{z}k_{\perp}^{2}$ where $k_{\perp}$ and $z$ are
transverse momentum and longitudinal momentum fraction carried by the daughter
gluon respectively. The inclusion of the kinematic constraint leads to a
modified GLR equation,

$\displaystyle\frac{\partial N(\eta,k_{\perp})}{\partial\eta}$
$\displaystyle=$
$\displaystyle\frac{\bar{\alpha}_{s}}{\pi}\int\frac{\mathrm{d}^{2}l_{\perp}}{l_{\perp}^{2}}N\left(\eta+\ln\left[\frac{k_{\perp}^{2}}{k_{\perp}^{2}+l_{\perp}^{2}}\right],l_{\perp}+k_{\perp}\right)-\frac{\bar{\alpha}_{s}}{\pi}\int_{0}^{k_{\perp}}\frac{\mathrm{d}^{2}l_{\perp}}{l_{\perp}^{2}}N(\eta,k_{\perp})-\bar{\alpha}_{s}N^{2}(\eta,k_{\perp}),$
(1)

with $\bar{\alpha}_{s}=\alpha_{s}N_{c}/\pi$, $\eta=\ln(x_{0}/x)$ and
$x_{0}=0.01$. The function $N(\eta,k_{\perp})$ is related to the normal TMD
gluon distribution $G(\eta,k_{\perp})$ through
$N(\eta,k_{\perp})=\frac{2\alpha_{s}\pi^{3}}{N_{c}S_{\perp}}G(\eta,k_{\perp})$
with $S_{\perp}$ being the transverse area of nucleon/nucleus. Converting the
above equation to the folded form of the GLR equation, it reads,

$\frac{\partial}{\partial\eta}\frac{N(x,k_{\perp})}{\Delta_{ns}(\eta,k_{\perp})}=\frac{\bar{\alpha}_{s}}{\pi}\int_{\Lambda_{\rm
cut}}\frac{\mathrm{d}^{2}l_{\perp}}{l_{\perp}^{2}}\frac{N\left(\eta+\ln\left[\frac{k_{\perp}^{2}}{k_{\perp}^{2}+l_{\perp}^{2}}\right],l_{\perp}+k_{\perp}\right)}{\Delta_{ns}(\eta,k_{\perp})}.$
(2)

where $\Delta_{ns}(\eta,k_{\perp})$ represents the probability of evolving
from $\eta_{0}$ to $\eta$ without resolvable branching. It is given by,

$\displaystyle\Delta_{ns}(\eta,k_{\perp})=\exp\left\\{-\bar{\alpha}_{s}\int^{\eta}_{\eta_{0}}d\eta^{\prime}\left[\ln\frac{k_{\perp}^{2}}{\Lambda_{\rm
cut}^{2}}+N(\eta^{\prime},k_{\perp})\right]\right\\},$ (3)

where the infrared cut off $\Lambda_{\rm cut}$ is the matter of choice about
what we classify as a resolvable emission. Emitted gluons with transverse
momentum $l_{\perp}<\Lambda_{\rm cut}$ are considered as the unresolvable
ones. And their contribution has been combined with the virtual correction to
cancel the infrared divergence. The resolvable branchings are defined as
emissions above this range. All order contributions from the virtual
correction and the unresolvable real emission are resummed into
$\Delta_{ns}(\eta,k_{\perp})$ which reduces to the non-Sudakov form factor
Kwiecinski  _et al._ (1996) in the dilute limit by neglecting the saturation
term. Eq. 2 can be converted into an integral form,

$N(\eta,k_{\perp})=N(\eta_{0},k_{\perp})\Delta_{ns}(\eta,k_{\perp})+\frac{\bar{\alpha}_{s}}{\pi}\int^{\eta}_{\eta_{0}}\mathrm{d}\eta^{\prime}\frac{\Delta_{ns}(\eta,k_{\perp})}{\Delta_{ns}(\eta^{\prime},k_{\perp})}\int_{\Lambda_{\rm
cut}}\frac{\mathrm{d}^{2}l_{\perp}}{l_{\perp}^{2}}N(\eta^{\prime}+\ln\left[\frac{k_{\perp}^{2}}{k_{\perp}^{2}+l_{\perp}^{2}}\right],l_{\perp}+k_{\perp}).$
(4)

It is evident that the kinematic constrained small $x$ equation is no longer a
local equation. Namely, the increase of gluon number density at rapidity
$\eta$ is driven by the gluon distribution at rapidity
$\eta+\ln\left[\frac{k_{\perp}^{2}}{k_{\perp}^{2}+l_{\perp}^{2}}\right]$
rather than that at the same rapidity $\eta$. The corresponding weighting
factor needs to be modified dramatically for the non-local case as shown
below.

###  II.1 Forward evolution

With these derived folded evolution equations, we are now ready to introduce
the Monte Carlo algorithm starting with the forward evolution case. For a
given initial condition $N(\eta_{i},k_{\perp,i})$, the first quantity to be
generated by the algorithm is the value of $\eta_{i+1}$. As it has been done
in (Shi  _et al._ , 2023), this task can be achieved by solving the equation,

$\mathcal{R}=\exp\left[-\bar{\alpha}_{s}\int^{\eta_{i+1}}_{\eta_{i}}\mathrm{d}\eta^{\prime}\left(\ln\frac{k_{\perp,i}^{2}}{\Lambda_{\rm
cut}^{2}}+N(\eta^{\prime},k_{\perp,i})\right)\right],$ (5)

where $\mathcal{R}$ is a random number distributed uniformly in the interval
[0,1]. Throughout this paper, we always use $\cal R$ to denote such a random
number. $N(\eta^{\prime},k_{\perp,i})$ is pre-generated by numerically solving
the GLR equation with the kinematic constraint.

In contrast to the DGLAP evolution, the unitarity is not preserved during the
course of the small $x$ evolution. The number of gluons increases after each
step of parton branching. The generated cascade thus needs to be re-weighted.
For instance, if one neglects the saturation effect and kinematic constraint
effect, the number of gluons which vanish due to the virtual correction and
the unresolved branching is proportional to
$\bar{\alpha}_{s}\int_{\Lambda_{\rm
cut}}^{k_{\perp,i}}\frac{\mathrm{d}l_{\perp}^{2}}{l_{\perp}^{2}}$, while the
number of gluons produced via the real correction is proportional to
$\bar{\alpha}_{s}\int_{\Lambda_{\rm
cut}}^{P_{\perp}}\frac{\mathrm{d}^{2}l_{\perp}}{l_{\perp}^{2}}$ where
$P_{\perp}$ is the UV cutoff, in the same rapidity interval. The weighting
function is given by the ratio of these two contributions ${\cal
W}(k_{\perp,i})=\ln(\frac{P_{\perp}^{2}}{\Lambda_{\rm
cut}^{2}})/\ln(\frac{k_{\perp,i}^{2}}{\Lambda_{\rm cut}^{2}})$.

Figure 1: Comparison of the gluon $k_{\perp}$ distributions reconstructed from
the forward evolution approach with the numerical solutions of the kinematic
constrained GLR equation at different rapidities. The results obtained with
the frozen boundary prescription and the fixed boundary prescription are shown
in the left and right plots respectively (Color online).

It is quite non-trivial to work out the correct weighting factor when the
kinematic constraint is implemented in the parton branching algorithm. Let us
first discuss the derivation of the weighting factor for the case of the fixed
boundary prescription. To work out the correct weighting coefficient, we first
write down the expression for the fraction of gluons at
$[\eta_{i+1},\eta_{i+1}+\delta\eta]$ that come form the branching between
$\eta_{i+1}$ and $\eta_{i}$,

$\displaystyle\delta\eta\frac{\partial}{\partial\eta_{i+1}}\left[\frac{\bar{\alpha}_{s}}{\pi}\int_{\eta_{i}}^{\eta_{i+1}}\mathrm{d}\eta^{\prime}\int_{\Lambda_{\rm
cut}}\frac{\mathrm{d}^{2}l_{\perp}}{l_{\perp}^{2}}e^{-\bar{\alpha}_{s}\int_{\eta_{i}}^{\eta^{\prime}}\mathrm{d}\eta\left[\ln\frac{k_{\perp,i}^{2}}{\Lambda_{\rm
cut}^{2}}+N(\eta,k_{\perp,i})\right]}\theta\left(\frac{1-z^{\prime}}{z^{\prime}}(k_{\perp,i}-l_{\perp})^{2}-l_{\perp}^{2}\right)\right]$
$\displaystyle=\delta\eta\frac{\bar{\alpha}_{s}}{\pi}\int^{{\rm
min}\left[P_{\perp},\sqrt{(k_{\perp,i}-l_{\perp})^{2}\frac{1-z}{z}}\right]}_{\Lambda_{\rm
cut}}\frac{\mathrm{d}^{2}l_{\perp}}{l_{\perp}^{2}}e^{-\bar{\alpha}_{s}\int_{\eta_{i}}^{\eta_{i+1}+\ln\frac{(k_{\perp,i}-l_{\perp})^{2}}{(k_{\perp,i}-l_{\perp})^{2}+l_{\perp}^{2}}}\mathrm{d}\eta\left[\ln\frac{k_{\perp,i}^{2}}{\Lambda_{\rm
cut}^{2}}+N(\eta,k_{\perp,i})\right]},$ (6)

with $z^{\prime}=x_{i+1}/x^{\prime}=\exp[\eta^{\prime}-\eta_{i+1}]$. The
kinematic constraint is imposed by the $\theta$-function. Note that the term
originating from the derivative acting on the integral boundary is equal to 0.
The entire contribution comes from the derivative acting on the
$\theta$-function. Meanwhile the fraction of gluons that leave from the
rapidity interval [$\eta_{i+1},\eta_{i+1}+\delta\eta$] due to the virtual
correction is,

$\displaystyle\delta\eta\frac{\partial
e^{-\bar{\alpha}_{s}\int_{\eta_{i}}^{\eta_{i+1}}\mathrm{d}\eta\left[\ln\frac{k_{\perp,i}^{2}}{\Lambda_{\rm
cut}^{2}}+N(\eta,k_{\perp,i})\right]}}{\partial\eta_{i+1}}=-\delta\eta\bar{\alpha}_{s}\left[\ln\frac{k_{\perp,i}^{2}}{\Lambda_{\rm
cut}^{2}}+N(\eta_{i+1},k_{\perp,i})\right]e^{-\bar{\alpha}_{s}\int_{\eta_{i}}^{\eta_{i+1}}\mathrm{d}\eta\left[\ln\frac{k_{\perp,i}^{2}}{\Lambda_{\rm
cut}^{2}}+N(\eta,k_{\perp,i})\right]}.$ (7)

For the non-local small $x$ evolution, one also needs the input for gluon
distribution beyond the small $x$ boundary $x_{0}=0.01$. There are two common
choices for the boundary conditions: i) the fixed boundary prescription,
$N(\eta<0,k_{\perp})=0$; ii) the frozen boundary prescription,
$N(\eta<0,k_{\perp})=N(\eta=0,k_{\perp})$. The weighting functions are thus
different for different rapidity boundary prescriptions.

For the fixed boundary prescription, the re-weighting function is given by

$\displaystyle{\cal
W}_{kc,1}(\eta_{i},\eta_{i+1};k_{\perp,i})=\frac{(\eta_{i+1}-\eta_{i})\int^{{\rm
min}\left[P_{\perp},\sqrt{\frac{1-z}{z}(k_{\perp,i}-l_{\perp})^{2}}\right]}_{\Lambda_{\rm
cut}}\frac{\mathrm{d}^{2}l_{\perp}}{l_{\perp}^{2}}e^{-\bar{\alpha}_{s}\int_{\eta_{i+1}}^{\eta_{i+1}+\ln\frac{(k_{\perp,i}-l_{\perp})^{2}}{(k_{\perp,i}-l_{\perp})^{2}+l_{\perp}^{2}}}\mathrm{d}\eta\left[\ln\frac{k_{\perp,i}^{2}}{\Lambda_{\rm
cut}^{2}}+N(\eta,k_{\perp,i})\right]}}{(\eta_{i+1}-\eta_{i})\ln\frac{k_{\perp,i}^{2}}{\Lambda_{\rm
cut}^{2}}+\int^{\eta_{i+1}}_{\eta_{i}}d\eta N(\eta,k_{\perp,i})}.$ (8)

Here, the values of $|l_{\perp}|$ and $\phi_{l}$ can be generated by solving
the following equation

$\displaystyle{\cal R}=\frac{1}{\cal
C}\frac{\bar{\alpha}_{s}}{\pi}\int_{\Lambda_{\rm
cut}}^{l_{\perp}}\frac{\mathrm{d}^{2}l^{\prime}_{\perp}}{l_{\perp}^{\prime
2}}\exp\left\\{-\bar{\alpha}_{s}\int_{\eta_{i}}^{\eta_{i+1}+\ln\frac{(k_{\perp,i}-l_{\perp}^{\prime})^{2}}{(k_{\perp,i}-l_{\perp}^{\prime})^{2}+l_{\perp}^{\prime
2}}}\mathrm{d}\eta\left[\ln\frac{k_{\perp,i}^{2}}{\Lambda_{\rm
cut}^{2}}+N(\eta,k_{\perp,i})\right]\right\\},$ (9) $\displaystyle{\cal
C}=\frac{\bar{\alpha}_{s}}{\pi}\int_{\Lambda_{\rm cut}}^{{\rm
min}[P_{\perp},\sqrt{(k_{\perp,i}-l_{\perp}^{\prime})^{2}\frac{1-z}{z}}]}\frac{\mathrm{d}^{2}l_{\perp}^{\prime}}{l_{\perp}^{\prime
2}}\exp\left\\{-\bar{\alpha}_{s}\int_{\eta_{i}}^{\eta_{i+1}+\ln\frac{(k_{\perp,i}-l_{\perp}^{\prime})^{2}}{(k_{\perp,i}-l_{\perp}^{\prime})^{2}+l_{\perp}^{\prime
2}}}\mathrm{d}\eta\left[\ln\frac{k_{\perp,i}^{2}}{\Lambda_{\rm
cut}^{2}}+N(\eta,k_{\perp,i})\right]\right\\},$ (10)

where ${\cal R}$ again is a random number and ${\cal C}$ is the normalization
factor ensuring that the r.h.s. of Eq. 9 resides in the region of $[0,1]$. In
the practical Monte Carlo implementation, a veto algorithm is used to be more
efficient. Once $|l_{\perp}|$ and $\phi_{l}$ are generated, $l$ and
$k_{\perp,i+1}$ then can be reconstructed subsequently. We repeat the
procedure outlined above until $\eta_{i+1}$ reach a minimal cut-off value
$\eta_{\rm min}$. Once the whole cascade is generated, we are able to
reconstruct the gluon $k_{\perp}$ distribution at arbitrary rapidity.

For the frozen boundary case, the weighting factor has to be modified to

$\mathcal{W}_{kc,2}(\eta_{i},\eta_{i+1};k_{\perp,i},k_{\perp,i+1})=\frac{(\eta_{i+1}-\eta_{i})\ln\frac{P_{\perp}^{2}}{\Lambda_{\rm
cut}^{2}}}{(\eta_{i+1}-\eta_{i})\ln\frac{k_{\perp,i}^{2}}{\Lambda_{\rm
cut}^{2}}+\int^{\eta_{i+1}}_{\eta_{i}}d\eta
N(\eta,k_{\perp,i})}\frac{{N(\eta_{i}+\ln\left[\frac{k_{\perp,i+1}^{2}}{k_{\perp,i+1}^{2}+l_{\perp}^{2}}\right],k_{\perp,i})}}{N(\eta_{i},k_{\perp,i})},$
(11)

and the radiated gluon transverse momentum $l_{\perp}$ is sampled solving the
following equation

${\cal R}=\frac{1}{\cal
C}\frac{\bar{\alpha}_{s}}{\pi}\int^{l_{\perp}}_{\Lambda_{\rm
cut}}\frac{\mathrm{d}^{2}l_{\perp}^{\prime}}{l_{\perp}^{\prime 2}},$ (12)

where the normalization factor for this case is given by ${\cal
C}=\frac{\bar{\alpha}_{s}}{\pi}\int^{P_{\perp}}_{\Lambda_{\rm
cut}}\frac{\mathrm{d}^{2}l_{\perp}^{\prime}}{l_{\perp}^{\prime 2}}$. The
$k_{\perp}$ distribution of the exchanged gluons that directly attaches to the
hard part can be reconstructed from the forward evolution algorithm described
above.

Using the recipes described above, we are now ready to generate parton
cascade. Following the conventional choice, we use the MV model McLerran and
Venugopalan (1994a, b) result as the initial condition at rapidity
$\eta_{0}=0$. Since we are interested in simulating events such as di-jet
production in eA collisions, it is suitable to utilize the Weisz$\ddot{a}$ke-
Williams (WW) gluon distribution as the initial condition Dominguez  _et al._
(2011). It is given by

$\displaystyle
N(\eta_{0},k_{\perp})=\int\frac{d^{2}r_{\perp}}{2\pi}e^{-ik_{\perp}\cdot
r_{\perp}}\frac{1}{r_{\perp}^{2}}\left(1-\exp\bigl{[}-\frac{1}{4}Q_{s0}^{2}r_{\perp}^{2}\ln(e+\frac{1}{\Lambda
r_{\perp}})\bigr{]}\right),$ (13)

with $Q_{s0}^{2}=1$ GeV2 and $\Lambda=0.24$ GeV. We explored the behavior of
the parton cascade with the both fixed boundary prescription and frozen
boundary prescription. From Fig. 1, one can see that the $k_{\perp}$
distribution obtained from the forward approach is in perfect agreement with
the numerical solutions of the kinematic constrained GLR equation for both
boundary conditions.

###  II.2 Backward evolution

We now turn to discuss how to implement the kinematic constraint in the
backward evolution which is far more efficient in generating initial state
parton shower as compared to the forward approach. The rapidity $\eta_{i+1}$
of gluon participating hard scattering is fixed by external kinematics.
$k_{\perp,i+1}$ at the rapidity $\eta_{i+1}$ can be sampled with the
distribution $N(\eta_{i+1},k_{\perp,i+1})$, which has to be determined
beforehand by numerically solving the evolution equation. The next step is to
generate $\eta_{i}$ using a modified non-Sudakov form factor.

The modified non-Sudakov form factor, $\Pi_{ns}$, can be related to the
forward non-Sudakov form factor $\Delta_{ns}$ and the gluon distribution $N$
as

$\displaystyle\Pi_{ns}(\eta_{i+1},\eta_{i};k_{\perp,i+1})=\frac{\Delta_{ns}(\eta_{i+1},k_{\perp,i+1})N(\eta_{i},k_{\perp,i+1})}{\Delta_{ns}(\eta_{i},k_{\perp,i+1})N(\eta_{i+1},k_{\perp,i+1})},$
(14)

which looks similar to that derived in our previous work Shi  _et al._ (2023).
However, one has to keep in mind that the gluon distributions appearing in the
above formula are obtained by solving the GLR equation with the kinematic
constraint.

On the other hand, the non-Sudakov factor can also be expressed as Shi  _et
al._ (2023),

$\displaystyle\Pi_{ns}(\eta_{i+1},\eta_{i};k_{\perp,i+1})=\exp\\!\\!\left[-\frac{\bar{\alpha}_{s}}{\pi}\int_{\eta_{i}}^{\eta_{i+1}}\\!\
\mathrm{d}\eta\int^{P_{\perp}}_{\Lambda_{\rm
cut}}\frac{\mathrm{d}^{2}l_{\perp}}{l_{\perp}^{2}}\frac{N\left(\eta\\!+\\!\ln\left[\frac{k_{\perp,i+1}^{2}}{k_{\perp,i+1}^{2}+l_{\perp}^{2}}\right]\\!,k_{\perp,i+1}+l_{\perp}\right)}{N(\eta,k_{\perp,i+1})}\right].$
(15)

Both non-Sudakov form factors can be equally well used to generate $\eta_{i}$
for a given $\eta_{i+1}$ by solving the following equation,

${\cal R}=\Pi_{ns}(\eta_{i+1},\eta_{i};k_{\perp,i+1}).$ (16)

The transverse momentum of the radiated gluon $l_{\perp}$ can be generated
according to

$\displaystyle{\cal R}=\frac{1}{\cal
C}\frac{\bar{\alpha}_{s}}{\pi}\int_{\Lambda_{\rm
cut}}^{l_{\perp}}\frac{\mathrm{d}^{2}l_{\perp}^{\prime}}{l_{\perp}^{\prime
2}}N\left(\eta_{i+1}+\ln\left[\frac{k_{\perp,i+1}^{2}}{k_{\perp,i+1}^{2}+l_{\perp}^{\prime
2}}\right],k_{\perp,i+1}+l_{\perp}^{\prime}\right),$ (17) $\displaystyle{\cal
C}=\frac{\bar{\alpha}_{s}}{\pi}\int_{\Lambda_{\rm
cut}}^{P_{\perp}}\frac{\mathrm{d}^{2}l_{\perp}^{\prime}}{l_{\perp}^{\prime
2}}N\left(\eta_{i+1}+\ln\left[\frac{k_{\perp,i+1}^{2}}{k_{\perp,i+1}^{2}+l_{\perp}^{\prime
2}}\right],k_{\perp,i+1}+l_{\perp}^{\prime}\right).$ (18)

Once again, $\cal R$ is a random number, $\cal C$ is the normalization factor
and a veto algorithm is employed in our practical implementation to make this
sampling procedure more efficient. Similar to the forward evolution case, the
generated event has to be re-weighted after each branching in the backward
evolution method as well. It is important to notice that the GLR equation with
the kinematic constraint is a non-local evolution equation when deriving the
weighting factor. The weighting factor associated with backward evolution is
the ratio of the fraction of gluons that appear from branching at the rapidity
$\eta_{i}+\ln\frac{k_{\perp,i+1}^{2}}{k_{\perp,i+1}^{2}+l_{\perp}^{2}}$ and
the fraction of gluons that vanish at the rapidity $\eta_{i}$ due to the
virtual correction and the fusion process. It reads,

$\mathcal{W}_{kc,\text{back}}(\eta_{i+1},\eta_{i};k_{\perp,i+1})=\frac{(\eta_{i+1}-\eta_{i})\ln\frac{k_{\perp,i}^{2}}{\Lambda_{\rm
cut}^{2}}+\int^{\eta_{i+1}}_{\eta_{i}}d\eta
N(\eta,k_{\perp,i})}{(\eta_{i+1}-\eta_{i})\ln\frac{P_{\perp}^{2}}{\Lambda_{\rm
cut}^{2}}}\frac{N(\eta_{i},k_{\perp,i})}{{N(\eta_{i}+\ln\left[\frac{k_{\perp,i+1}^{2}}{k_{\perp,i+1}^{2}+l_{\perp}^{2}}\right],k_{\perp,i})}}.$
(19)

The procedure outlined above is repeated until $\eta_{i}$ is smaller than
$\eta_{0}$. The last step of the simulation is to construct four momenta of
the radiated gluons. Note that the minus component of the $t$-channel gluon’s
four momentum can only be reconstructed after the full cascade has been
generated. By going from the last $t$-channel gluon (closest to the nucleus),
which has the vanishing minus component, forward in the cascade to the hard
scattering process, the true minus component of the $t$-channel gluons are
constructed. In Fig. 2, we compare gluon $k_{\perp}$ distribution at different
rapidities generated from backward evolution to the numerical solutions of the
GLR equation with the kinematic constraint. The perfect match between gluon
$k_{\perp}$ distributions obtained from the backward approach and by
numerically solving the kinematic constrained GLR has been found.

Figure 2: Comparison of the gluon $k_{\perp}$ distributions obtained from the
backward approach with the numerical solutions of the GLR equation at
different rapidities (Color online).

##  III  $k_{t}$ resummation in the small $x$ limit

Our ultimate goal is to build a parton shower generator for simulating events
in eA collisions at EIC. The hard scattering processes occurring in eA
collisions often involve multiple scales. For instance, loosely speaking,
there are three well separated scales in the back-to-back di-jet production:
the center mass of energy $\sqrt{s}$, the invariant mass of the di-jet $Q$,
and the total transverse momentum of the di-jet system $k_{\perp}$. To improve
the convergence of the pertubative series, the two type large logarithms
$\alpha_{s}\ln\left({s}/{Q^{2}}\right)$ and
$\alpha_{s}\ln^{2}\left({Q^{2}}/{k_{\perp}^{2}}\right)$ arise in the high
order calculations of the di-jet production cross section have to be summed to
all orders. The summation of the logarithm contribution
$\alpha_{s}\ln\left({s}/{Q^{2}}\right)$ is achieved by solving the small $x$
evolution equation, while the logarithm contribution
$\alpha_{s}\ln^{2}\left({Q^{2}}/{k_{\perp}^{2}}\right)$ can be resummed by
means of the CS equation. A unified framework that allows us to resum both
large logarithms simultaneously in a consistent way have been developed in a
sequence of papers Zhou (2016); Xiao  _et al._ (2017); Zhou (2019). The
evolved small $x$ gluon TMD can be expressed as the convolution of the Sudakov
form factor and the renormalized dipole amplitudes. It has been stressed in
Refs. Xiao  _et al._ (2017); Zhou (2019) that at small $x$, gluon TMDs only
can be matched onto dipole scattering amplitudes rather than the normal gluon
PDFs in the collinear factorization. We notice that such a joint resummation
formalism has been studied in the various different context Mueller  _et al._
(2013a, b); Zheng  _et al._ (2014); Stasto  _et al._ (2018); Marquet  _et al._
(2020); van Hameren  _et al._ (2019, 2021); Shi  _et al._ (2022); Zhao  _et
al._ (2021); Boussarie  _et al._ (2021); Hentschinski (2021); Caucal  _et al._
(2021); Taels  _et al._ (2022); Caucal  _et al._ (2022); Goda  _et al._
(2022); Benić  _et al._ (2022); Al-Mashad  _et al._ (2022); Goda (2023);
Caucal  _et al._ (2023); Ganguli  _et al._ (2023); Goda  _et al._ (2023); van
Hameren  _et al._ (2023).

To simulate hard scattering processes involving multiple scales in a parton
shower generator, it is necessary to develop a Monte Carlo branching algorithm
to effectively resum both types of logarithms through an iteration procedure.
The essential observation that enables the computer implementation of the
joint resummation is described as the following. In the backward approach, the
evolution starts from the final $t$-channel gluon with the most negative
virtual mass-squared, which participates in the hard process. As a parton
cascade develops towards the backward direction, the virtual mass of the
$t$-channel gluon decreases by radiating soft gluons with the longitudinal
momentum fraction $1-z\rightarrow 0$. This first stage of the evolution is
described by the CS equation and the renormalization group equation which
resum the double leading $k_{t}$ logarithm and the single leading $k_{t}$
logarithm respectively. When the virtual mass of the $t$-channel gluon goes
down to the scale which is of the order of saturation scale, we should perform
the small $x$ evolution. The precise value of this scale should be fixed by
fitting the output of the cascade to the experimental data. During the course
of the small $x$ evolution, the virtual mass of the $t$-channel gluon stops
monotonously decreasing, whereas its longitudinal momentum fraction increases
rapidly until the small $x$ evolution initial boundary is reached. In this
second stage of the evolution, the development of parton cascade is mainly
driven by the radiated gluons that carry the large longitudinal momentum
fraction $1-z\rightarrow 1$. Therefore, the Monte Carlo algorithm based on the
GLR equation should be applied to generate the parton branching at this stage.

To simulate the first stage of the evolution, our primary task is to derive a
folded version of the CS equation and the renormalization group equation. To
this end, we write down the CS equation in the momentum space,

$\displaystyle\frac{\partial
N(\mu^{2},\zeta^{2},\eta,k_{\perp})}{\partial\ln\zeta^{2}}=\frac{\bar{\alpha}_{s}}{2\pi}\int^{\zeta}_{0}\frac{d^{2}l_{\perp}}{l_{\perp}^{2}}\left[N(\mu^{2},\zeta^{2},\eta,k_{\perp}+l_{\perp})-N(\mu^{2},\zeta^{2},\eta,k_{\perp})\right]\
.$ (20)

which can be converted into the conventional expression of the CS equation
Collins  after making the Fourier transform up to the leading logarithm
accuracy. Here, $\mu$ is the factorization scale, and $\zeta$ is a scale
introduced to regularize the light cone divergence. The factorization scale
dependence of the gluon TMD in the saturation regime is described by the
normal renormalization group equation Xiao  _et al._ (2017),

$\displaystyle\frac{\partial
N(\mu^{2},\zeta^{2},\eta,k_{\perp})}{\partial\ln\mu^{2}}=\bar{\alpha}_{s}\left[\beta_{0}-\frac{1}{2}\ln\frac{\zeta^{2}}{\mu^{2}}\right]N(\mu^{2},\zeta^{2},\eta,k_{\perp})\
.$ (21)

with $\beta_{0}=\frac{11}{12}-\frac{N_{f}}{6N_{c}}$ and $N_{f}=3$ in this
work. By choosing the factorization scale $\mu$ to be $\zeta$, one can combine
the CS equation and the renormalization group equation together. The combined
evolution equation reads,

$\displaystyle\frac{\partial N(Q^{2},\eta,k_{\perp})}{\partial\ln
Q^{2}}=\frac{\bar{\alpha}_{s}}{2\pi}\int^{Q}_{0}\frac{d^{2}l_{\perp}}{l_{\perp}^{2}}\left[N(Q^{2},\eta,k_{\perp}+l_{\perp})-N(Q^{2},\eta,k_{\perp})\right]+\bar{\alpha}_{s}\beta_{0}N(Q^{2},\eta,k_{\perp}),$
(22)

where $N(Q^{2},\eta,k_{\perp})\equiv
N(\mu^{2}=Q^{2},\zeta^{2}=Q^{2},\eta,k_{\perp})$. Following the standard
procedure, the above evolution equation can be cast into a folded equation,

$\displaystyle\frac{\partial}{\partial\ln
Q^{2}}\frac{N(Q^{2},\eta,k_{\perp})}{\Delta_{s}(Q^{2})}=\frac{\bar{\alpha}_{s}}{2\pi}\int_{\Lambda_{\rm
cut}}^{Q}\frac{d^{2}l_{\perp}}{l_{\perp}^{2}}\frac{N(Q^{2},\eta,k_{\perp}+l_{\perp})}{\Delta_{s}(Q^{2})},$
(23)

with the Sudakov form factor being given by,

$\displaystyle{\Delta_{s}(Q^{2})}=\exp\left[-\int_{Q_{0}^{2}}^{Q^{2}}\frac{dt}{t}\frac{\bar{\alpha}_{s}(t)}{2}\left(\ln\frac{t}{\Lambda_{\rm
cut}^{2}}-2\beta_{0}\right)\right].$ (24)

The Sudakov form factor is simply the probability of evolving from $Q_{0}$ to
$Q$ without branching. Eq. 23 can be integrated to give an integral equation
for $N(Q^{2},\eta,k_{\perp})$ in terms of the gluon TMD at the initial scale
$Q_{0}$:

$\displaystyle
N(Q^{2},\eta,k_{\perp})=N(Q_{0}^{2},\eta,k_{\perp}){\Delta_{s}(Q^{2})}+\int^{{Q}^{2}}_{Q^{2}_{0}}\frac{dt}{t}\frac{\Delta_{s}(Q^{2})}{\Delta_{s}(t)}\frac{\bar{\alpha}_{s}(t)}{2\pi}\int_{\Lambda_{\rm
cut}}^{Q}\frac{d^{2}l_{\perp}}{l_{\perp}^{2}}N(t,\eta,k_{\perp}+l_{\perp}).$
(25)

With the derived folded CS and renormalization group equation, we are ready to
introduce the Monte Carlo implementation of the $k_{t}$ resummation formulated
in the framework of the CGC effective theory.

###  III.1 Forward evolution

To have a consistency check, we first present the formulation of the forward
evolution scheme. The combined CS and renormalization group equation can be
solved using the forward evolution approach. We lay out the main procedures in
the following.

For a given virtuality scale $Q_{i}$, either after several steps of evolution
or at the initial condition, we first generate the value of a higher
virtuality scale $Q_{i+1}$, where the next branching occurs. Following the
conventional method, this can be achieved by solving the following equation,

$\displaystyle\mathcal{R}=\exp\left[-\int_{Q_{i}^{2}}^{Q_{i+1}^{2}}\frac{dt}{t}\bar{\alpha}_{s}(t)\left(\frac{1}{2}\ln\frac{t}{\Lambda_{\rm
cut}^{2}}-\beta_{0}\right)\right].$ (26)

where the argument of the running coupling $\alpha_{s}$ is simply chosen to be
the virtual mass squared.

Once $Q_{i+1}$ is generated, the transverse momentum of the radiated gluon,
$l_{\perp,i+1}$, can be determined according to the following equation

$\displaystyle\mathcal{R}=\frac{1}{\cal C}\int_{\Lambda_{\rm
cut}}^{l_{\perp,i+1}}\frac{d^{2}{l^{\prime}}_{\perp}}{{l^{\prime}}_{\perp}^{2}},$
(27)

where the normalization factor reads ${\cal C}=\int_{\Lambda_{\rm
cut}}^{Q_{i+1}}\frac{d^{2}{l^{\prime}}_{\perp}}{{l^{\prime}}_{\perp}^{2}}$.
The four momenta of the radiated gluon and the $t$-channel gluon can be
determined from the momentum conservation and the on-shell condition. We will
discuss the reconstruction of kinematics in more details in the next
subsection.

The generated cascade needs to be re-weighted. This is because that the
unitary is no longer preserved beyond the leading double logarithm
approximation. We have included the leading single logarithm contribution in
the algorithm employed here, which leads to the increase of gluon number
density after each splitting. The weighting factor is given by,

${\cal W}_{\rm
CS}(Q^{2}_{i+1},Q^{2}_{i})=\frac{\int^{Q_{i+1}^{2}}_{Q_{i}^{2}}\frac{dt}{t}\alpha_{s}(t)\ln\frac{t}{\Lambda_{\rm
cut}^{2}}}{\int^{Q_{i+1}^{2}}_{Q_{i}^{2}}\frac{dt}{t}\alpha_{s}(t)\left[\ln\frac{t}{\Lambda_{\rm
cut}^{2}}-2\beta_{0}\right]}.$ (28)

If the single logarithm contribution associated with the $\beta_{0}$ term in
the denominator is neglected, the weighting factor reduces to 1. With these
re-weighted parton cascades, one can reconstruct the $t$-channel gluon
$k_{\perp}$ distribution at different scales and compare with the analytical
and numerical solutions of Eq. 22.

It is straightforward to numerically solve Eq. 22, while the analytical
solution of Eq. 22 can also be easily obtained in the impact parameter space.
After Fourier transforming back to the momentum space, the evolved gloun TMD
distribution reads,

$N(Q^{2},\eta,k_{\perp})=\int\frac{d^{2}b_{\perp}}{(2\pi)^{2}}e^{ik_{\perp}\cdot
b_{\perp}}e^{-S(\mu_{b}^{2},Q^{2})}\int d^{2}l_{\perp}e^{-il_{\perp}\cdot
b_{\perp}}N(\eta,l_{\perp}),$ (29)

where $N(\eta,l_{\perp})$ is the gluon distribution evolved with the GLR
equation, or the initial condition computed in the MV model. The Sudakov
factor at one loop level in the impact parameter ($b_{\perp}$) space consists
of a perturbative part and a non-perturbative part. It is given by

$\displaystyle
S(\mu_{b}^{2},Q^{2})=S_{pert}(\mu_{b*}^{2},Q^{2})+S_{NP}(b_{\perp}^{2},Q^{2}).$
(30)

The perturbative Sudakov factor reads

$\displaystyle
S_{pert}(\mu_{b*}^{2},Q^{2})=\frac{N_{c}}{2\pi}\int^{Q^{2}}_{\mu_{b*}^{2}}\frac{d\mu^{2}}{\mu^{2}}\alpha_{s}(\mu)\left[\ln\frac{Q^{2}}{\mu^{2}}-2\beta_{0}\right],$
(31)

where $\mu_{b*}^{2}$ is defined as
$\mu_{b*}^{2}=4e^{-2\gamma_{E}}/b_{\perp*}^{2}$, with
$b_{\perp*}=\frac{b_{\perp}}{\sqrt{1+b_{\perp}^{2}/b^{2}_{\max}}}$ and
$b_{\max}=1.5\text{ GeV}^{-1}$. To compare with the Monte Carlo result on the
same footing, we simply neglect the non-perturbative Sudakov factor $S_{NP}$
in the numerical calculation. The behaviour at large $b_{\perp}$ is regulated
by $N(\eta,b_{\perp})$ which is the Fourier transform of $N(\eta,l_{\perp})$.
In this work, we use the one-loop running coupling which reads

$\alpha_{s}(\mu^{2})=\frac{1}{\beta_{0}\frac{N_{c}}{\pi}\ln(\mu^{2}/\Lambda^{2}_{\rm
QCD})},$ (32)

with $\Lambda^{2}_{\rm QCD}=0.0578$ GeV2.

We present the $t$-channel gluon $k_{\perp}$ distribution constructed from the
generated parton cascade and compare it with the numerical solution of the CS-
renormalization group equation for the fixed coupling case in the left panel
of Fig. 3. In our estimation, the MV model is employed to provide with the
gluon distribution at the initial scale $Q_{0}=$3 GeV. In the formulation of
TMD evolution, all soft-radiated gluons carry exactly zero longitudinal
momentum fraction. In contrast, all radiated soft gluons carry finite
longitudinal momentum fraction in the parton branching algorithm. This
presents an important advantage of the Monte Carlo method comparing with the
conventional analytical approach. Keeping longitudinal momentum conservation
exactly in parton splitting process is often crucial to correctly account for
phenomenology near the threshold region Liu  _et al._ (2022). However, to make
the comparisons in a consistent way, we didn’t change the longitudinal
momentum fraction of the $t$-channel gluon after each branching in our
algorithm. In the right panel of Fig. 3, we compare the Monte Carlo simulation
result with both the numerical solution of the CS-renormalization equation and
the analytical solution for the running coupling case at the scale $Q=13$ GeV.
It is clear to see from the right panel of Fig. 3 that our algorithm yields
the same $k_{\perp}$ distribution as the numerical result. On the other hand,
it differs from the analytical approach result. Such discrepancy is expected
because the non-perturbative part of the CS kernel is treated differently in
the analytical approach. In addition, the argument of the running coupling
used in the parton branching algorithm and the numerical solution is the hard
scale $Q$, whereas the scale of running coupling is $\mu_{b}$ in the
analytical approach. Since the analytical result can describe the relevant
phenomenology very well, one should use it as guidance to model the non-
perturbative part of the Sudakov factor which will be introduced in the Monte
Carlo algorithm in the future work. Alternatively, one could also use a
relatively large infrared cutoff value $\Lambda_{\rm cut}$ to mimic the effect
of the non-perturbative Sudakov factor. We leave this for a future study.

Figure 3: Comparison of the gluon $k_{\perp}$ distributions obtained from the
forward evolution approach with the numerical solutions of the combined CS-
renormalization group equation at different scales(Color online). The left
panel: the fixed coupling case. The right panel: the running coupling
case(Color online).

###  III.2 Backward evolution

In this subsection, we will outline the essential steps of Monte Carlo
implementation for the backward evolution based on the folded CS-
renormalization group evolution equation. Unlike the forward evolution which
can be considered as a way of solving the evolution equation, the evolved
parton distributions have to be pre-generated and are used to guide the
backward evolution. In the most parton branching algorithm, the $k_{t}$
resummation is achieved by using the modified Sudakov factor incorporating the
collinear Parton Distributions Functions (PDFs). However, in the saturation
regime, the $k_{t}$ resummation has to be formulated in terms of the
unintegrated gluon distribution. The main procedures are summarized as
follows.

The modified Sudakov factor in the backward evolution approach is different
from that in the forward evolution approach. It reads

$\Pi_{s}(Q_{i+1},Q_{i};k_{\perp,i+1})=\frac{{\Delta_{s}(Q_{i+1}^{2})}N(Q^{2}_{i},\eta,k_{\perp,i+1})}{{\Delta_{s}(Q_{i}^{2})}N(Q^{2}_{i+1},\eta,k_{\perp,i+1})}.$
(33)

An alternative way to compute the modified Sudakov factor is given by

$\displaystyle\Pi_{s}(Q_{i+1},Q_{i};k_{\perp,i+1})$ $\displaystyle=$
$\displaystyle\exp\left[-\int_{Q_{i}^{2}}^{Q_{i+1}^{2}}\frac{dt}{t}\frac{\bar{\alpha}_{s}(t)}{2\pi}\int^{\sqrt{t}}_{\Lambda_{\rm
cut}}\frac{d^{2}l_{\perp}}{l_{\perp}^{2}}\frac{N(t,\eta,k_{\perp,i+1}+l_{\perp})}{N(t,\eta,k_{\perp,i+1})}\right].$
(34)

It describes the probability for gluon evolving backward from $Q_{i+1}$ to
$Q_{i}$ without branching. The transverse momentum dependent gluon
distribution appearing in Eq. 33 and Eq. 34 has to be pre-generated by
numerically solving the combined CS-renormalization group equation.

The backward evolution starts from the $t$-channel gluon with the highest
virtuality $Q_{i}$. The hard scale of the partonic scattering process is
denoted as $Q_{i+1}$. We first have to sample $k_{\perp,i+1}$ according to the
following distribution

$\mathcal{R}=\frac{1}{\cal C}\int^{k_{\perp,i+1}}_{\Lambda_{\rm
cut}}d^{2}k^{\prime}_{\perp}N(Q_{i+1}^{2},\eta,k_{\perp}^{\prime}),$ (35)

with ${\cal C}=\int^{Q_{i+1}}_{\Lambda_{\rm
cut}}d^{2}k^{\prime}_{\perp}N(Q_{i+1}^{2},\eta,k_{\perp}^{\prime})$ being the
normalization factor. The rapidity $\eta$ is fixed by external kinematics. The
next quantity to be generated by the parton cascade algorithm is the value of
virtuality $Q_{i}$.

Following the standard backward evolution strategy, $Q_{i}$ is obtained using
the backward type Sudakov factor. We can sample a $Q_{i}$ by solving the
following equation,

$\displaystyle{\cal R}=\Pi_{s}(Q_{i+1},Q_{i};k_{\perp,i+1}).$ (36)

As the virtual mass of $(i+1)$th $t$-channel gluon, $Q_{i}$ also serves as the
hard probe scale at which the $i$th $t$-channel gluon’s transverse momentum is
measured. The transverse momentum of the radiated gluon $l_{\perp,i}$ is thus
sampled solving the following equation

$\displaystyle\mathcal{R}=\frac{1}{\cal C}\int^{l_{\perp,i}}_{\Lambda_{\rm
cut}}\frac{d^{2}l^{\prime}_{\perp}}{l^{\prime
2}_{\perp}}N(Q_{i}^{2},\eta,k_{\perp,i+1}+l^{\prime}_{\perp}),$ (37)
$\displaystyle{\cal C}=\int^{Q_{i}}_{\Lambda_{\rm
cut}}\frac{d^{2}l^{\prime}_{\perp}}{l^{\prime
2}_{\perp}}N(Q_{i}^{2},\eta,k_{\perp,i+1}+l^{\prime}_{\perp}).$ (38)

The longitudinal momentum fraction of the radiated gluon is determined through
the on-shell condition,

$|Q_{i}^{2}|\approx\frac{z_{i}l_{\perp,i}^{2}}{1-z_{i}}+|k_{\perp,i+1}^{2}|,$
(39)

which is valid in the strong ordering region
$|Q_{i-1}^{2}|\ll|Q_{i}^{2}|\ll|Q_{i+1}^{2}|$. The minus component of the
emitted gluon can be fixed accordingly. The $i$th $t$-channel gluon’s
transverse momentum is trivially obtained:
$k_{\perp,i}=k_{\perp,i+1}-l_{\perp,i}$. The virtual mass $Q_{i-1}$ of the
$i$th $t$-channel gluon is computed with Eq. 36. However, $t$-channel gluons’
four momenta can be determined only after the whole cascade is generated. The
minus component of the $t$-channel gluon that is directly attached to nuclear
target is set to be 0. From this initial condition, the four momenta of
$t$-channel gluons are retrospectively reconstructed by momentum conservation.

As argued in the previous subsection, the generated event has to be re-
weighted after each branching since the unitary is not preserved in the single
leading logarithm accuracy level. In the backward evolution approach, the re-
weighting function reads

${\cal W}_{\rm
CS,back}(Q^{2}_{i+1},Q^{2}_{i})=\frac{\int^{Q_{i+1}^{2}}_{Q_{i}^{2}}\frac{dt}{t}\alpha_{s}(t)\left[\ln\frac{t}{\Lambda_{\rm
cut}^{2}}-2\beta_{0}\right]}{\int^{Q_{i+1}^{2}}_{Q_{i}^{2}}\frac{dt}{t}\alpha_{s}(t)\ln\frac{t}{\Lambda_{\rm
cut}^{2}}}\ .$ (40)

We repeat the procedure outlined above until $Q_{i}^{2}$ reach a minimal cut-
off scale at which TMD evolution stops. The TMD evolution is driven by the
soft gluon radiations which carry the vanishing longitudinal momentum fraction
$1-z_{i}\rightarrow 0$. In the practical Monte Carlo implementation, the cut-
off is chosen to be $|Q_{i}^{2}|>|l_{\perp,i}^{2}|+|k_{\perp,i+1}^{2}|$, or
equivalently $z_{i}>0.5$. Meanwhile, $|Q_{i}^{2}|$ is also required to be
larger than the satuartion scale $Q_{s}^{2}$. If these two conditions can not
be met simultaneously, we terminate the TMD evolution, and start the backward
small $x$ evolution.

We test the backward evolution algorithm against the numerical method as shown
in Fig. 4. The MV model result is applied at the initial scale $Q_{0}$=3 GeV.
The gluon $k_{\perp}$ distribution at high scale $Q=13$ GeV is obtained by
numerically solving the combined CS-renormalization group equation. The
cascade is generated starting from the scale $Q=13$ GeV and evolve down to the
initial scale with the backward approach. The $t$-channel gluon $k_{\perp}$
distribution reconstructed from the cascade is compared with the numerical
results at different scales. Gluon $k_{\perp}$ distributions are presented in
the left panel of Fig. 4 for the fixed coupling case, and in the right panel
of Fig. 4 for the running coupling case. It is evident that the $k_{\perp}$
distributions obtained from the Monte Carlo method is the same as the
numerical results. We conclude that the backward evolution algorithm pass this
consistency check as expected.

Figure 4: Comparison of the gluon $k_{\perp}$ distributions obtained from the
backward approach with the numerical solutions of the CS-renormalzation group
equation at different scales (Color online).

##  IV Conclusion

In this work, we extended the small $x$ initial state parton branching
algorithm developed in the previous paper to include the kinematic constraint
effect. In the small $x$ limit, the kinematic constraint leads to stronger
suppression of soft gluon emissions than that caused by the angular ordering
along the chain. The coherent branching effect is thus effectively implemented
in the parton branching algorithm once the kinematic constraint is imposed.
This is a nontrivial extension in the sense that the weighting factor and the
way of sampling radiated gluon’s transverse momenta are drastically altered.
The $t$-channel gluon $k_{\perp}$ distributions constructed from both the
forward scheme and the backward scheme are shown to reproduce the numerical
solutions of the kinematic constrained GLR equation.

We also formulated a parton branching algorithm that enables us to resum large
$k_{t}$ logarithms at small $x$ logarithms following a two-step evolution
picture. The cascade first develops by radiating soft gluons that carry
vanishing longitudinal momentum fractions in the backward approach
description. At this first stage of the evolution, the parton branching is
simulated with the Sudakov factor which we obtained from the folded CS
equation and the renormalization group equation. The transverse momentum-
dependent gluon distribution instead of gluon PDF is used to guide the
evolution path toward the most populated regions of $(Q^{2},k_{\perp})$. When
the virtual mass of the $t$-channel gluon is dominated by its transverse
momentum or is of the order of saturation scale, the parton branching starts
being generated according to the non-Sudakov form factor derived from the
small $x$ evolution equation. The joint $k_{t}$ and small $x$ resummation thus
has been achieved in the Monte Carlo simulation by implementing such two-step
evolution. Our study represents an important step towards practical
applications of the parton shower generator in simulating scattering processes
that involve multiple well-separated hard scales, such as di-jet production in
eA collisions at EIC. The next step is to construct a full hadron-level Monte
Carlo generator with the hadronization being performed using multi-purpose
generators such as PYTHIA Sjostrand  _et al._ (2006). We also plan to
integrate the algorithm into eHIJING framework Ke  _et al._ (2023) aiming at
the simulation of events in eA collisions for the whole $x$ range accessible
at EIC in the future.

Acknowledgments: We thank Hai-tao Li and Shan-shan Cao for helpful
discussions. This work has been supported by the National Natural Science
Foundation of China under Grant No. 1217511. Y.S. is supported by the China
Postdoctoral Science Foundation under Grant No. 2022M720082. S.Y.W. is also
supported by the Taishan fellowship of Shandong Province for junior
scientists.

## References

  * Accardi  _et al._ (2016) A. Accardi _et al._ , Eur. Phys. J. A 52, 268 (2016), arXiv:1212.1701 [nucl-ex] .
  * Abdul Khalek  _et al._ (2022) R. Abdul Khalek _et al._ , Nucl. Phys. A 1026, 122447 (2022), arXiv:2103.05419 [physics.ins-det] .
  * Ellis  _et al._ (2011) R. K. Ellis, W. J. Stirling, and B. R. Webber,  _QCD and collider physics_, Vol. 8 (Cambridge University Press, 2011).
  * Buckley  _et al._ (2011) A. Buckley _et al._ , Phys. Rept. 504, 145 (2011), arXiv:1101.2599 [hep-ph] .
  * Höche (2015) S. Höche, in  _Theoretical Advanced Study Institute in Elementary Particle Physics: Journeys Through the Precision Frontier: Amplitudes for Colliders_ (2015) pp. 235–295, arXiv:1411.4085 [hep-ph] .
  * Byer  _et al._ (2022) D. Byer, V. Khachatryan, H. Gao, I. Akushevich, A. Ilyichev, C. Peng, A. Prokudin, S. Srednyak, and Z. Zhao, (2022), arXiv:2210.03785 [hep-ph] .
  * Altarelli and Parisi (1977) G. Altarelli and G. Parisi, Nucl. Phys. B 126, 298 (1977).
  * Marchesini and Webber (1991) G. Marchesini and B. R. Webber, Nucl. Phys. B 349, 617 (1991).
  * Marchesini and Webber (1992) G. Marchesini and B. R. Webber, Nucl. Phys. B 386, 215 (1992).
  * Jung and Salam (2001) H. Jung and G. P. Salam, Eur. Phys. J. C 19, 351 (2001), arXiv:hep-ph/0012143 .
  * Jung  _et al._ (2010a) H. Jung _et al._ , Eur. Phys. J. C 70, 1237 (2010a), arXiv:1008.0152 [hep-ph] .
  * Andersen and Smillie (2010a) J. R. Andersen and J. M. Smillie, JHEP 01, 039 (2010a), arXiv:0908.2786 [hep-ph] .
  * Andersen and Smillie (2010b) J. R. Andersen and J. M. Smillie, Phys. Rev. D 81, 114021 (2010b), arXiv:0910.5113 [hep-ph] .
  * Jung  _et al._ (2010b) H. Jung _et al._ (CASCADE), Eur. Phys. J. C 70, 1237 (2010b), arXiv:1008.0152 [hep-ph] .
  * Andersen and Smillie (2011) J. R. Andersen and J. M. Smillie, JHEP 06, 010 (2011), arXiv:1101.5394 [hep-ph] .
  * Lipatov  _et al._ (2020) A. V. Lipatov, M. A. Malyshev, and S. P. Baranov, Eur. Phys. J. C 80, 330 (2020), arXiv:1912.04204 [hep-ph] .
  * van Hameren (2018) A. van Hameren, Comput. Phys. Commun. 224, 371 (2018), arXiv:1611.00680 [hep-ph] .
  * Hautmann  _et al._ (2022) F. Hautmann, M. Hentschinski, L. Keersmaekers, A. Kusina, K. Kutak, and A. Lelek, Phys. Lett. B 833, 137276 (2022), arXiv:2205.15873 [hep-ph] .
  * Kuraev  _et al._ (1977) E. A. Kuraev, L. N. Lipatov, and V. S. Fadin, Sov. Phys. JETP 45, 199 (1977).
  * Balitsky and Lipatov (1978) I. I. Balitsky and L. N. Lipatov, Sov. J. Nucl. Phys. 28, 822 (1978).
  * Ciafaloni (1988) M. Ciafaloni, Nucl. Phys. B 296, 49 (1988).
  * Catani  _et al._ (1990a) S. Catani, F. Fiorani, and G. Marchesini, Nucl. Phys. B 336, 18 (1990a).
  * Catani  _et al._ (1990b) S. Catani, F. Fiorani, and G. Marchesini, Phys. Lett. B 234, 339 (1990b).
  * Marchesini (1995) G. Marchesini, Nucl. Phys. B 445, 49 (1995), arXiv:hep-ph/9412327 .
  * Lipatov  _et al._ (2023) A. V. Lipatov, G. I. Lykasov, and M. A. Malyshev, Phys. Rev. D 107, 014022 (2023), arXiv:2211.03727 [hep-ph] .
  * Lipatov and Malyshev (2023) A. V. Lipatov and M. A. Malyshev, (2023), arXiv:2305.04005 [hep-ph] .
  * Shi  _et al._ (2023) Y. Shi, S.-Y. Wei, and J. Zhou, Phys. Rev. D 107, 016017 (2023), arXiv:2211.07174 [hep-ph] .
  * Gribov  _et al._ (1983) L. V. Gribov, E. M. Levin, and M. G. Ryskin, Phys. Rept. 100, 1 (1983).
  * Balitsky (1996) I. Balitsky, Nucl. Phys. B 463, 99 (1996), arXiv:hep-ph/9509348 .
  * Kovchegov (1999) Y. V. Kovchegov, Phys. Rev. D 60, 034008 (1999), arXiv:hep-ph/9901281 .
  * Jalilian-Marian  _et al._ (1997) J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, Nucl. Phys. B 504, 415 (1997), arXiv:hep-ph/9701284 .
  * Jalilian-Marian  _et al._ (1998) J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, Phys. Rev. D 59, 014014 (1998), arXiv:hep-ph/9706377 .
  * Iancu  _et al._ (2001) E. Iancu, A. Leonidov, and L. D. McLerran, Nucl. Phys. A 692, 583 (2001), arXiv:hep-ph/0011241 .
  * Ferreiro  _et al._ (2002) E. Ferreiro, E. Iancu, A. Leonidov, and L. McLerran, Nucl. Phys. A 703, 489 (2002), arXiv:hep-ph/0109115 .
  * Kwiecinski  _et al._ (1996) J. Kwiecinski, A. D. Martin, and P. J. Sutton, Z. Phys. C 71, 585 (1996), arXiv:hep-ph/9602320 .
  * Kwiecinski  _et al._ (1997) J. Kwiecinski, A. D. Martin, and A. M. Stasto, Phys. Rev. D 56, 3991 (1997), arXiv:hep-ph/9703445 .
  * Deak  _et al._ (2019) M. Deak, K. Kutak, W. Li, and A. M. Staśto, Eur. Phys. J. C 79, 647 (2019), arXiv:1906.09062 [hep-ph] .
  * Zhou (2016) J. Zhou, JHEP 06, 151 (2016), arXiv:1603.07426 [hep-ph] .
  * Xiao  _et al._ (2017) B.-W. Xiao, F. Yuan, and J. Zhou, Nucl. Phys. B 921, 104 (2017), arXiv:1703.06163 [hep-ph] .
  * Zhou (2019) J. Zhou, Phys. Rev. D 99, 054026 (2019), arXiv:1807.00506 [hep-ph] .
  * Liu  _et al._ (2022) H.-y. Liu, X.-h. Liu, Y. Shi, D.-x. Zheng, and J. Zhou, Phys. Rev. D 106, 036026 (2022), arXiv:2204.00262 [hep-ph] .
  * Avsar  _et al._ (2011) E. Avsar, A. M. Stasto, D. N. Triantafyllopoulos, and D. Zaslavsky, JHEP 10, 138 (2011), arXiv:1107.1252 [hep-ph] .
  * Iancu  _et al._ (2015a) E. Iancu, J. D. Madrigal, A. H. Mueller, G. Soyez, and D. N. Triantafyllopoulos, Phys. Lett. B 744, 293 (2015a), arXiv:1502.05642 [hep-ph] .
  * Iancu  _et al._ (2015b) E. Iancu, J. D. Madrigal, A. H. Mueller, G. Soyez, and D. N. Triantafyllopoulos, Phys. Lett. B 750, 643 (2015b), arXiv:1507.03651 [hep-ph] .
  * Lappi and Mäntysaari (2016) T. Lappi and H. Mäntysaari, Phys. Rev. D 93, 094004 (2016), arXiv:1601.06598 [hep-ph] .
  * Ducloué  _et al._ (2019) B. Ducloué, E. Iancu, A. H. Mueller, G. Soyez, and D. N. Triantafyllopoulos, JHEP 04, 081 (2019), arXiv:1902.06637 [hep-ph] .
  * Ducloué  _et al._ (2020) B. Ducloué, E. Iancu, G. Soyez, and D. N. Triantafyllopoulos, Phys. Lett. B 803, 135305 (2020), arXiv:1912.09196 [hep-ph] .
  * Zheng and Zhou (2019) D.-X. Zheng and J. Zhou, JHEP 11, 177 (2019), arXiv:1906.06825 [hep-ph] .
  * McLerran and Venugopalan (1994a) L. D. McLerran and R. Venugopalan, Phys. Rev. D 49, 2233 (1994a), arXiv:hep-ph/9309289 .
  * McLerran and Venugopalan (1994b) L. D. McLerran and R. Venugopalan, Phys. Rev. D 49, 3352 (1994b), arXiv:hep-ph/9311205 .
  * Dominguez  _et al._ (2011) F. Dominguez, C. Marquet, B.-W. Xiao, and F. Yuan, Phys. Rev. D 83, 105005 (2011), arXiv:1101.0715 [hep-ph] .
  * Mueller  _et al._ (2013a) A. H. Mueller, B.-W. Xiao, and F. Yuan, Phys. Rev. Lett. 110, 082301 (2013a), arXiv:1210.5792 [hep-ph] .
  * Mueller  _et al._ (2013b) A. H. Mueller, B.-W. Xiao, and F. Yuan, Phys. Rev. D 88, 114010 (2013b), arXiv:1308.2993 [hep-ph] .
  * Zheng  _et al._ (2014) L. Zheng, E. C. Aschenauer, J. H. Lee, and B.-W. Xiao, Phys. Rev. D 89, 074037 (2014), arXiv:1403.2413 [hep-ph] .
  * Stasto  _et al._ (2018) A. Stasto, S.-Y. Wei, B.-W. Xiao, and F. Yuan, Phys. Lett. B 784, 301 (2018), arXiv:1805.05712 [hep-ph] .
  * Marquet  _et al._ (2020) C. Marquet, S.-Y. Wei, and B.-W. Xiao, Phys. Lett. B 802, 135253 (2020), arXiv:1909.08572 [hep-ph] .
  * van Hameren  _et al._ (2019) A. van Hameren, P. Kotko, K. Kutak, and S. Sapeta, Phys. Lett. B 795, 511 (2019), arXiv:1903.01361 [hep-ph] .
  * van Hameren  _et al._ (2021) A. van Hameren, P. Kotko, K. Kutak, and S. Sapeta, Phys. Lett. B 814, 136078 (2021), arXiv:2010.13066 [hep-ph] .
  * Shi  _et al._ (2022) Y. Shi, L. Wang, S.-Y. Wei, and B.-W. Xiao, Phys. Rev. Lett. 128, 202302 (2022), arXiv:2112.06975 [hep-ph] .
  * Zhao  _et al._ (2021) Y.-Y. Zhao, M.-M. Xu, L.-Z. Chen, D.-H. Zhang, and Y.-F. Wu, Phys. Rev. D 104, 114032 (2021), arXiv:2105.08818 [hep-ph] .
  * Boussarie  _et al._ (2021) R. Boussarie, H. Mäntysaari, F. Salazar, and B. Schenke, JHEP 09, 178 (2021), arXiv:2106.11301 [hep-ph] .
  * Hentschinski (2021) M. Hentschinski, Phys. Rev. D 104, 054014 (2021), arXiv:2107.06203 [hep-ph] .
  * Caucal  _et al._ (2021) P. Caucal, F. Salazar, and R. Venugopalan, JHEP 11, 222 (2021), arXiv:2108.06347 [hep-ph] .
  * Taels  _et al._ (2022) P. Taels, T. Altinoluk, G. Beuf, and C. Marquet, JHEP 10, 184 (2022), arXiv:2204.11650 [hep-ph] .
  * Caucal  _et al._ (2022) P. Caucal, F. Salazar, B. Schenke, and R. Venugopalan, (2022), arXiv:2208.13872 [hep-ph] .
  * Goda  _et al._ (2022) T. Goda, K. Kutak, and S. Sapeta, (2022), arXiv:2210.16084 [hep-ph] .
  * Benić  _et al._ (2022) S. Benić, O. Garcia-Montero, and A. Perkov, Phys. Rev. D 105, 114052 (2022), arXiv:2203.01685 [hep-ph] .
  * Al-Mashad  _et al._ (2022) M. A. Al-Mashad, A. van Hameren, H. Kakkad, P. Kotko, K. Kutak, P. van Mechelen, and S. Sapeta, JHEP 12, 131 (2022), arXiv:2210.06613 [hep-ph] .
  * Goda (2023) T. Goda, Acta Phys. Polon. Supp. 16, 18 (2023).
  * Caucal  _et al._ (2023) P. Caucal, F. Salazar, B. Schenke, T. Stebel, and R. Venugopalan, (2023), arXiv:2304.03304 [hep-ph] .
  * Ganguli  _et al._ (2023) I. Ganguli, A. van Hameren, P. Kotko, and K. Kutak, (2023), arXiv:2306.04706 [hep-ph] .
  * Goda  _et al._ (2023) T. Goda, K. Kutak, and S. Sapeta, (2023), arXiv:2305.14025 [hep-ph] .
  * van Hameren  _et al._ (2023) A. van Hameren, H. Kakkad, P. Kotko, K. Kutak, and S. Sapeta, (2023), arXiv:2306.17513 [hep-ph] .
  * (74) J. C. Collins, Foundations of perturbative QCD(Cambridge University Press, Cambridge, 2011) .
  * Sjostrand  _et al._ (2006) T. Sjostrand, S. Mrenna, and P. Z. Skands, JHEP 05, 026 (2006), arXiv:hep-ph/0603175 .
  * Ke  _et al._ (2023) W. Ke, Y.-Y. Zhang, H. Xing, and X.-N. Wang, (2023), arXiv:2304.10779 [hep-ph] .