  Multi-particle azimuthal cumulants, often used to study collective flow in
high-energy heavy-ion collisions, have recently been applied in small collision
systems such as $pp$ and $p$+A to extract the second-order azimuthal harmonic
flow $v_2$. Recent observation of four-, six- and eight-particle cumulants with
"correct sign" $c_2\{4\}<0, c_2\{6\}>0, c_2\{8\}<0$ and approximate equality of
the inferred single-particle harmonic flow, $v_2\{4\}\approx v_2\{6\}\approx
v_2\{8\}$, have been used as strong evidence for a collective emission of all
soft particles produced in the collisions. We show that these relations in
principle could be violated due to the non-Gaussianity in the event-by-event
fluctuation of flow and/or non-flow. Furthermore, we show, using $pp$ events
generated with the PYTHIA model, that $c_2\{2k\}$ obtained with standard
cumulant method are dominated by non-flow from dijets. An alternative cumulant
method based on two or more $\eta$-separated subevents is proposed to suppress
the dijet contribution. The new method is shown to be able to recover a flow
signal as low as 4\% imposed on the PYTHIA events, independently of how the
event activity class is defined. Therefore the subevent cumulant method offers
a more robust way of studying collectivity based on the existence of long-range
azimuthal correlations between multiple distinct $\eta$ ranges. The prospect of
using the subevent cumulants to study collective flow in A+A collisions, in
particular its longitudinal dynamics, is discussed.
