Let $(M^n,g)$ be a Riemannian manifold that admit a unit Killing vector field $X$. i.e., $\mathscr{L}_Xg=0$. Is it possible that there exist a smooth function $f$ on $M$ such that $X=\mathrm{grad}f$?
A good reference will be very appreciated.
It is possible for this to happen (for example, the unit translation vector fields in $\mathbb{R}^n$), but it does not necessarily hold (for example, there are unit Killing vector fields on the stanard Riemannian $3$-sphere, but, obviously, none of them could be gradients, since they never vanish and the $3$-sphere is compact).
I'm not sure what fact or theory you want a reference for.
The Killing condition $\nabla_a X^\flat_b + \nabla_b X^\flat_a = 0$, together with the gradient condition (really, the closedness condition for a $1$-form) $\nabla_a X^\flat_b - \nabla_b X^\flat_a = 0$, imply that $X$ is covariantly constant, $\nabla_a X^b = 0$. Note that I'm using the relation $X^\flat(-) = \langle X, - \rangle$ to define the $1$-form $X^\flat$.
The existence of covariantly constant vectors implies various restrictive conditions on the geometry of the Riemannian manifold $(M^n, g)$. If there's any particular restriction that you are interested in, you may want to refine your question.