In this article, we extend previous studies based on CMB spherical harmonic phases (SHP) to examine the validity of the hypothesis that the temperature field of the CMB is consistent with a Gaussian random field (GRF). The null hypothesis is that the corresponding CMB SHP are independent and identically distributed in terms of a uniform distribution in the interval [0, 2$\pi$] \citep{1986ApJ...304...15B,2013rossmanith}. We devise a new model-independent method where we use ordered and non-parametric Rao's statistic, based on sample arc-lengths to comprehensively test uniformity and independence of SHP. We performed our analysis on the scales limited by spherical harmonic modes $\le$ 128, to restrict ourselves to signal-dominated regions. To find the non-uniform or dependent sets of SHP, we calculate the statistic for the data and 10000 Monte Carlo simulated uniformly random sets of SHP and use 0.05 and 0.001 $\alpha$ levels to distinguish between statistically significant and highly significant detections. We first establish the performance of our method using simulated Gaussian, non-Gaussian CMB temperature maps, along with observed non-Gaussian 100 and 143 GHz Planck channel maps. We find that our method performs efficiently and accurately in detecting phase correlations generated in all of the non-Gaussian simulations and observed foreground contaminated 100 and 143 GHz Planck channel temperature maps. We apply our method on the Planck satellite mission's final released CMB temperature anisotropy maps- COMMANDER, SMICA, NILC, and SEVEM along with WMAP 9 year released ILC map. We report that SHP corresponding to some of the $m$-modes is non-uniform, some of the $\ell$ mode SHP and neighboring mode pair SHP are correlated in cleaned CMB maps. The detection of non-uniformity or correlation in the SHP indicates the presence of non-Gaussian signals in the foreground minimized CMB maps., Comment: 42 pages, 23 figures, 6 tables