X-ray pulsar navigation(XPNAV) is an attractive method for autonomous navigation of deep space in the future. Currently, techniques for estimating the phase of X-ray pulsar radiation involve the maximization of the general non-convex object functions based on the average profile from the epoch folding method. This results in the suppression of useful information and highly complex computation. In this paper, a new maximum likelihood(ML) phase estimation method that directly utilizes the measured time of arrivals(TOAs) is presented. The X-ray pulsar radiation will be treated as a cyclo-stationary process and the TOAs of the photons in a period will be redefined as a new process, whose probability distribution function is the normalized standard profile of the pulsar. We demonstrate that the new process is equivalent to the generally used Poisson model. Then, the phase estimation problem is recast as a cyclic shift parameter estimation under the ML estimation, and we also put forward a parallel ML estimation method to improve the ML solution. Numerical simulation results show that the estimator described here presents a higher precision and reduces the computational complexity compared with currently used estimators. X-ray pulsar navigation (XPNAV) is an attractive method for autonomous navigation of deep space in the future. Currently, techniques for estimating the phase of X-ray pulsar radiation involve the maximization of the general non-convex object functions based on the average profile from the epoch folding method. This result in the suppression of useful information and highly complex computation. In this paper, a new maximum likelihood (ML) phase estimation method that directly takes the measured time of arrivals (TOAs) is presented. -ray pulsar radiation will be treated as a cyclo-stationary process and the TOAs of the photons in a period will be redefined as a new process, whose probability distribution function is the normalized standard profile of the pulsar. We demonstrate that the new process is Then, the phase estimation problem is recast as a cyclic shift parameter estimation under the ML estimation, and we also put forward ap Numerical simulation results show that the estimator described here presents a higher precision and reduces the computational complexity compared with currently used estimators.