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\left(p_T\right)$ and $v_3\left(p_T\right)$. 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\text{fm}/\text{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}\left(p_T\right)$ 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.