  Correlations of two flow harmonics $v_n$ and $v_m$ via three- and
four-particle cumulants are measured in 13 TeV $pp$, 5.02 TeV $p$+Pb, and 2.76
TeV peripheral Pb+Pb collisions with the ATLAS detector at the LHC. The goal is
to understand the multi-particle nature of the long-range collective phenomenon
in these collision systems. The large non-flow background from dijet production
present in the standard cumulant method is suppressed using a method of
subevent cumulants involving two, three and four subevents separated in
pseudorapidity. The results show a negative correlation between $v_2$ and $v_3$
and a positive correlation between $v_2$ and $v_4$ for all collision systems
and over the full multiplicity range. However, the magnitudes of the
correlations are found to depend strongly on the event multiplicity, the choice
of transverse momentum range and collision system. The relative correlation
strength, obtained by normalisation of the cumulants with the $\langle
v_n^2\rangle$ from a two-particle correlation analysis, is similar in the three
collision systems and depends weakly on the event multiplicity and transverse
momentum. These results based on the subevent methods provide strong evidence
of a similar long-range multi-particle collectivity in $pp$, $p$+Pb and
peripheral Pb+Pb collisions.
