  In this paper, we investigate the hydrodynamic collectivity in
high-multiplicity events of proton-proton collisions at $\sqrt{s}=$ 13 TeV,
using iEBE-VISHNU hybrid model with three different initial conditions, namely,
HIJING, super-MC and TRENTo. With properly tuned parameters, hydrodynamic
simulations with each initial model give reasonable descriptions of the
measured two-particle correlations, including the integrated and $p_{\rm
T}$-differential flow for all charged and identified hadrons. However, the
hydrodynamic simulations fail to describe the negative value of the
four-particle cumulant $c_2^v\{4\}$ as measured in experiments. Further
investigations show that the non-linear response between the elliptic flow
$v_2$ and the initial eccentricity $\varepsilon_2$ becomes significant in the
small p-p systems. This leads to a large deviation from linear eccentricity
scaling and generates additional flow fluctuations, which results in a positive
$c_2^v\{4\}$ even with a negative $c_2^\varepsilon\{4\}$ from the initial
state. We also presented the first hydrodynamic calculations of multi-particle
mixed harmonic azimuthal correlations in p-p collisions, such as normalized
asymmetric cumulant $nac_n\{3\}$, normalized Symmetric-Cumulant,
$nsc_{2,3}\{4\}$ and $nsc_{2,4}\{4\}$. Although many qualitative features are
reproduced by the hydrodynamic simulations with chosen parameters, the measured
negative $nsc_{2,3}\{4\}$ cannot be reproduced. The failure of the description
of negative $c_2\{4\}$ and $nsc_{2,3}\{4\}$ triggers the question on whether
hydrodynamics with a fundamentally new initial state model could solve this
puzzle, or hydrodynamics itself might not be the appreciated mechanism of the
observed collectivity in p-p collisions at the LHC.
