  We numerically study for the first time the nonlinear GLR-MQ evolution
equations for nuclear parton distribution function (nPDFs) to next-to-leading
order accuracy and quantify the impact of gluon recombination at small $x$.
Using the nCTEQ15 nPDFs as input, we confirm the importance of the nonlinear
corrections for small $x \lesssim 10^{-3}$, whose magnitude increases with a
decrease of $x$ and an increase of the atomic number $A$. We find that at
$x=10^{-5}$ and for heavy nuclei, after the upward evolution from $Q_0=2$ GeV
to $Q=10$ GeV, the quark singlet $\Omega(x,Q^2)$ and the gluon $G(x,Q^2)$
distributions become reduced by $9-15$%, respectively. The relative effect is
much stronger for the downward evolution from $Q_0=10$ GeV to $Q=2$ GeV, where
we find that $\Omega(x,Q^2)$ is suppressed by 40%, while $G(x, Q^2)$ is
enhanced by 140%. These trends propagate into the $F_2^A(x,Q^2)$ nuclear
structure function and the $F_L^A(x,Q^2)$ longitudinal structure function,
which after the downward evolution become reduced by 45% and enhanced by 80%,
respectively. Our analysis indicates that the nonlinear effects are most
pronounced in $F_L^A(x,Q^2)$ and are already quite sizable at $x \sim 10^{-3}$
for heavy nuclei. We have checked that our conclusions very weakly depend on
the choice of input nPDFs. In particular, using the EPPS21 nPDFs as input, we
obtain quantitatively similar results.
