For all one dimensional bound states (be it relativistic or non-relativistic), user Gonenc has explained clearly why $\langle p \rangle =0$.
However, I want to know whether $\langle p \rangle=0$ is also true for 3-dimensional bound states. Are there any criteria such as relativistic or non-relativistic, symmetric or asymmetric, etc. under which $\langle p \rangle=0$ holds? Justification for your answer and reference to your answer would be much appreciated.
In general of course not (not even in $1d$). However, the following observations hold.
If the wavefunction is purely real, then $\langle p_i\rangle=0$ for $i=x,y,z$. This is easy since $p_i\mapsto -i\hbar d/dx_i$ so that (for instance) \begin{align} \langle p_x\rangle = -i\hbar \int dx\,dy\,dz\, \psi^*(x,y,z)\frac{\partial}{\partial x}\psi(x,y,z)\, . \end{align} If the right hand side is non-zero, then it must be necessarily purely imaginary since $\psi\in \mathbb{R}$ and $\partial \psi/\partial x\in \mathbb{R}$. However, the average value of a Hermitian operator is necessarily real; thus the hypothesis that the right hand side is non-zero is false. This holds for any Cartesian component in any dimension.
The argument fails for a number of important states: for instance a coherent state $\vert\alpha\rangle$ where $\alpha$ has an imaginary part. Note that a coherent state is not an eigenstate of the harmonic oscillator hamiltonian.
It is also easy to find linear combinations of harmonic oscillator states or particle-in-a-box states for which $\langle p\rangle\ne 0$. $\frac{1}{\sqrt{2}}e^{-i E_1t/\hbar}\vert 1\rangle +e^{-iE_2t/\hbar}\frac{i}{\sqrt{2}}\vert 2\rangle$ is such an example. Again such a state is not an eigenstate of $\hat H$ (either for h.o. of particle-in-a-box.)
Broadly speaking, since $\langle p\rangle \propto \frac{d}{dt}\langle x\rangle$, $\langle p\rangle$ will be $0$ whenever $\langle x\rangle$ does not depend on time. This is guaranteed to be so if $\psi$ is a solution to the time-independent Schrodinger equation since in this case $\vert\psi\vert^2$ will not depend on $t$ so that $\langle x\rangle = \int dx x\vert\psi\vert^2$ will also not depend on $t$.
The converse is not true: there are cases where, even if $\langle x\rangle$ is constant (and thus time-independent), $\vert \psi\vert^2$ may depend on $t$. For instance the combination of harmonic oscillator states $a_0\vert 0\rangle e^{-i E_0t/\hbar} + a_2\vert 2\rangle e^{-i E_2t/\hbar}$ with $\vert a_0\vert^2+\vert a_2\vert^2=1$ has $\langle x\rangle =0$ but the probability density $\vert \psi\vert^2$ is still time dependent. Note that for this case $\langle p\rangle=0$, irrespective of $a_0$ and $a_2$.
