  A detailed study of multi-particle azimuthal correlations is presented using
$pp$ data at $\sqrt{s}=5.02$ and 13 TeV, and $p$+Pb data at
$\sqrt{s_{\rm{NN}}}=5.02$ TeV, recorded with the ATLAS detector at the LHC. The
azimuthal correlations are probed using four-particle cumulants $c_{n}\{4\}$
and flow coefficients $v_n\{4\}=(-c_{n}\{4\})^{1/4}$ for $n=2$ and 3, with the
goal of extracting long-range multi-particle azimuthal correlation signals and
suppressing the short-range correlations. The values of $c_{n}\{4\}$ are
obtained as a function of the average number of charged particles per event,
$\left\langle N_{\rm{ch}} \right\rangle$, using the recently proposed
two-subevent and three-subevent cumulant methods, and compared with results
obtained with the standard cumulant method. The three-subevent method is found
to be least sensitive to short-range correlations, which originate mostly from
jets with a positive contribution to $c_{n}\{4\}$. The three-subevent method
gives a negative $c_{2}\{4\}$, and therefore a well-defined $v_2\{4\}$, nearly
independent of $\left\langle N_{\rm{ch}} \right\rangle$, which provides direct
evidence that the long-range multi-particle azimuthal correlations persist to
events with low multiplicity. Furthermore, $v_2\{4\}$ is found to be smaller
than the $v_2\{2\}$ measured using the two-particle correlation method, as
expected for long-range collective behavior. Finally, the measured values of
$v_2\{4\}$ and $v_2\{2\}$ are used to estimate the number of sources relevant
for the initial eccentricity in the collision geometry.
