  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.
