  We compute single and double inclusive gluon distributions in classical
Yang-Mills simulations of proton-lead collisions and extract the associated
transverse momentum dependent Fourier harmonics $v_2(p_T)$ and $v_3(p_T)$.
Gluons have a large $v_2$ in the initial state, while odd harmonics such as
$v_3$ vanish identically at the initial time $\tau=0^{+}$. By the time $\tau
\lesssim 0.4\,{\rm fm/c}$ final state effects in the classical Yang-Mills
evolution generate a non-zero $v_{3}$ and only mildly modify the gluon $v_{2}$.
Unlike hydrodynamic flow, these momentum space anisotropies are uncorrelated
with the global spatial anisotropy of the collision. A principal ingredient for
the generation of $v_2$ and $v_3$ in this framework is the event-by-event
breaking of rotational invariance in domains the size of the inverse of the
saturation scale $Q_s$. In contrast to our findings in p+Pb collisions
Yang-Mills simulations of lead-lead collisions generate much smaller values of
$v_{2,3} (p_T)$ and additional collective flow effects are needed to explain
experimental data. This is because the locally generated anisotropy due to the
breaking of rotational invariance is depleted with the increase in the number
of uncorrelated domains.
