  We revisit the problem of small Bjorken-$x$ evolution of the gluon and
flavor-singlet quark helicity distributions in the shock wave ($s$-channel)
formalism. Earlier works on the subject in the same framework resulted in an
evolution equation for the gluon field-strength $F^{12}$ and quark "axial
current" ${\bar\psi}\gamma^+\gamma^5\psi$ operators (sandwiched between the
appropriate light-cone Wilson lines) in the double-logarithmic approximation
(DLA: summing powers of $\alpha_s\,\ln^2(1/x)$ with $\alpha_s$ the strong
coupling constant). In this work, we observe that an important mixing of the
above operators with another gluon operator, $\overleftarrow{D}^i\,{D}^i$, also
sandwiched between the light-cone Wilson lines (with the repeated index $i=1,2$
summed over), was missing in the previous works. This operator has the physical
meaning of the sub-eikonal (covariant) phase: its contribution to helicity
evolution is shown to be proportional to another sub-eikonal operator,
$D^i-\overleftarrow{D}^i$, which is related to the Jaffe-Manohar polarized
gluon distribution. In this work we include this operator into small-$x$
helicity evolution, and construct a novel evolution mixing all three operators
($D^i-\overleftarrow{D}^i$, $F^{12}$, and ${\bar\psi}\gamma^+\gamma^5\psi$),
generalizing the previous results. We also construct closed DLA evolution
equations in the large-$N_c$ and large-$N_c\& N_f$ limits, with $N_c$ and $N_f$
the numbers of quark colors and flavors, respectively. Solving the large-$N_c$
equations numerically we obtain the following small-$x$ asymptotics of the
quark and gluon helicity distributions $\Delta \Sigma$ and $\Delta G$, along
with the $g_1$ structure function, \[\Delta\Sigma(x,Q^2)\sim\Delta G(x,Q^2)\sim
g_1(x,Q^2)\sim\left(\frac{1}{x}\right)^{3.66\,\sqrt{\frac{\alpha_s\,N_c}{2\pi}}},\]
in complete agreement with the earlier work by Bartels, Ermolaev and Ryskin.
