The Kuramoto model is universally applicable to many weakly interacting systems, such as Josephson junctions. However, in many cases couplings which include higher harmonics are needed, such as those relating to electro-chemical oscillators and phi-Josephson junctions. In this study, we consider a global mean field synchronization model, called the nonlinear mean field Kuramoto model, where a finite number of phase oscillators are coupled via a global mean field which is a nonlinear function of the Kuramoto mean field. As observed first by Komarov and Pikovsky [1], identical oscillators with initial phase drawn randomly from a uniform distribution form two asymmetrical clusters under such a model. The distribution obeys a scaling law of N^1/2, showing that such a symmetry breaking is related to initial fluctuations. We attempt to find the source in the asymmetry through analytical method as well as numerical method. Analytically, we conducted linear stability analysis, and reformulated the Watanabe-Strogatz dynamical equations in the complex plane. The singularity in the Moebius transformation is analytically located, which provides a candidate for the singularity of the dynamics that we seek. The numerical simulation confirms the location predicted by theory, but it also confirms that statistical tools can't be applied to study the location of the singularity based on only the initial condtions, because the singularity cannot be obtained without integration. During the investigation, further properties of Watanabe-Strogatz variables and of the nonlinear mean field Kuramoto model were discovered to indicate several future paths to take for uncovering this process of symmetry breaking.

[1] M. Komarov, A. Pikovsky, Finite-size-induced transitions to synchrony in oscillator ensembles with nonlinear global coupling, Phys. Rev. E, v. 92, 020901(R) 2015