  In this paper, we propose a framework to investigate the collective dynamics
in ensembles of globally coupled phase oscillators when higher-order modes
dominate the coupling. The spatiotemporal properties of the attractors in
various regions of parameter space are analyzed. Furthermore, a detailed linear
stability analysis proves that the stationary symmetric distribution is only
neutrally stable in the marginal regime which stems from the generalized
time-reversal symmetry. Moreover, the critical parameters of the transition
among various regimes are determined analytically by both the Ott-Antonsen
method and linear stability analysis, the transient dynamics are further
revealed in terms of the characteristic curves method. Finally, for the more
general initial condition the symmetric dynamics could be reduced to a rigorous
three-dimensional manifold which shows that the neutrally stable chaos could
also occur in this model for particular parameter. Our theoretical analysis and
numerical results are consistent with each other, which can help us understand
the dynamical properties in general system with higher-order harmonics
couplings.
