Phase synchronization of weakly coupled limit cycle oscillators are related to the stability of the zero solution of the reduced-order dynamics of phase differences, represented by a systems of differential equations on a hypertorus. Using Rantzer's density function, a dual form of Lyapunov function, we propose a method to certify almost global stability of an equilibrium on a hypertorus. We show that the proposed method can certify robustness of phase synchronization of all-to-all and weakly coupled limit cycle oscillators with respect to disturbances in phases. The method leverages sum of squares polynomial optimization to construct the certification function.