Distance function on complete Riemannian manifold.

Let $C_p$ be the cut locus of $p$. $C_p$ contains two type of points: (i) points $q$ such that there exist at least two minimizing geodesics from $p$ to $q$; (ii) points $q$ that are conjugate to $p$. It can be proven that $C_p$ is closed and also a null set (i.e. if $(h,U)$ is a chart, then $h (U \cap C_p)$ is a null set with respect to the usual Lebesgue measure).

There is an alternative description of $C_p$. Let $I_p$ be the connected component of $\{ v \in T_p M \mid \Bbb d _v \exp_p \text{ is an isomorphism} \}$ containing $0 \in T_p M$. Let $c_p = \partial I_p$, the boundary of $I_p$. Then $C_p = \exp _p (c_p)$.

Finally, to answer your question: it is well known that $d(p, \cdot)$ is smooth on $M \setminus (C_p \cup \{p\})$ (that $d(p, \cdot)$ cannot be smooth in $p$ is obvious: choose normal coordinates around $p$ and think of what happens to the Euclidean norm around $0$). The square $d(p, \cdot)^2$ is smooth on $M \setminus C_p$.

For further details on the cut locus, see for instance sub-chapter III.2 of "Riemannian Geometry" by Isaac Chavel (I have a feeling, though, that you will not enjoy his style). A clearer exposition can be found in sub-chapter VIII.7 of vol. 2 of "Foundations of Differential Geometry" by Kobayashi & Nomizu (conjugate points are introduced in VIII.3).

A brief exposition of these concepts, together with a proof of the differentiability of $d(p, \cdot)$ can be found in theorem 3.1 of the lecture notes of a course given by Zuoqin Wang at the University of Science and Technology of China.