Expectation values of $(x,y,z)$ in the $|n\ell m\rangle$ state of hydrogen? Expectation values of $(x,y,z)$ in the $| n\ell m\rangle$ state of hydrogen?
Does anyone know of a quick way of finding this (if there is even one)? Can I somehow use the relation that:
$$\langle r\rangle ~=~ \frac{a_0}{2}(3n^2 - \ell(\ell+1)),$$ 
or do I just have to brute force and use properties of Laguerre polynomials and spherical harmonics and what not?
 A: You're confusing the expectation values of the vector $\mathbf{r}=(x,y,z)$ and its magnitude $r=\sqrt{x^2+y^2+z^2}$. Because the hydrogen hamiltonian is parity invariant, all its eigenfunctions are chosen to have a definite parity. This means that the expectation value $\langle\mathbf{r}\rangle$ will always be zero because each component is the integral of an odd function times an even probability distribution function.
The second expectation value you mention,
$$\langle n,l,m|r|n,l,m\rangle=\int \mathrm{d}x \mathrm{d}y \mathrm{d}z \psi_{n,l,m}^\ast(x,y,z)\sqrt{x^2+y^2+z^2}\psi_{n,l,m}(x,y,z),$$
is quite different, since you're not counting the direction of the distances you add. This is the same as comparing the averages of $x$ and $|x|$ over some 1D probability distribution function $p(x)$. The average of $x$ will be zero if $p(x)$ is even, while the average of $|x|$ will only vanish is $p$ is a delta function, concentrated at the origin. Thus $\langle|x|\rangle$ can be quite large while $\langle x\rangle$ is zero because of cancellations.

This is not to say that expectation values of $\mathbf{r}$ are not interesting, but one must simply be more careful. The fact that the diagonal matrix elements vanish says that the eigenstates have no permanent dipole moment - which of course they can't as they are eigenstates of an isotropic system. What doesn't vanish, however, are the transition matrix elements,
$$\langle n',l',m'|\mathbf{r}|n,l,m\rangle,$$
which play an important role in the hydrogen atom's interaction with radiation. (Specifically, they control the leading-order interaction energy, which is the dipole coupling $-q\mathbf{r}\cdot\mathbf{E}$.) Because the components of $\mathbf{r}$ are themselves spherical harmonics of degree 1 (or linear combinations of them), the matrix element above will involve the spherical integral
$$\int\mathrm{d}\Omega \, Y_{l'm'}^\ast Y_{1,\mu} Y_{lm}$$
for $\mu=-1,0,1$. 
Because the spherical harmonics have a rich orthogonality structure, only a few of these integrals will survive. Specifically, you need $\Delta l=|l-l'|=1$ and $m'=m+\mu$, which has the physical content of the dipole selection rules: for a dipole coupling to radiation, only S-P, P-D, D-F, ... transitions are allowed. Further, only a few spatial orientations are allowed: $m$ can increase or decrease by one (which happens if $|\mu|=1$, and corresponds to circularly polarized light in either direction, with the electric field rotating in the $x,y$ plane) or stay the same (when $m=0$ and the light is linearly polarized along $z$).
A: I think if you guess, you'll get the right answer!
Edit: What I've written below is not quite right; it ignores the a priori possibility that two degenerate energy eigenstates could have different values of $k$, and a superposition of them might not satisfy $\langle \hat x\rangle = 0$.  So you do need a bit more information.  You can modify the argument below by re-instating $y$ and $z$, and considering the parity operation $\vec{x} \to -\vec{x}$.  Then use the property that $Y_{lm}(-\vec{x}) = (-1)^l Y_{lm}(\vec{x})$.

But you can also calculate it without knowing anything about spherical harmonics etc.  Let's focus on $x$, and ask what $\langle \psi | \hat x | \psi \rangle$ is for some (unspecified) energy eigenstate $\psi$.  First note that the Hamiltonian for the problem is invariant under $x \to -x$, so without loss of generality, we can assume that $\psi$ is an eigenstate of the corresponding operator, i.e., $\psi(-x) = (-1)^k \psi(x)$ for $k = 0$ or $1$.  Now let's just calculate (for simplicity of notation, I will suppress the integration over, and dependence on, $y$ and $z$; it makes no difference):
\begin{eqnarray*}
\langle\psi|\hat x|\psi\rangle &=& ~\int_{-\infty}^{\infty}dx\, \bar\psi(x)\,x\,\psi(x)\\[2ex]
&\stackrel{x~ :=~ -x'}=& \int_{-\infty}^{\infty}dx'(-1)^{2k} \bar\psi(x')(-x')\psi(x')\\[2ex]
&=&~ -\langle\psi|\hat x|\psi\rangle
\end{eqnarray*}
We conclude that the expectation value is $0$.  Obviously the same goes for $y$ and $z$.
(I think there is a subtlety in the above argument: it only holds for normalisable states, i.e., the bound states.  Unbound energy eigenstates can only be delta-function normalised, and the above breaks down.)
A: Yes. The vector $(x,y,z)$ transforms as a spin-1 tensor. So by the rules of the addition of angular momentum, its expectation value is only nonzero for $L=1$ states of the Hydrogen atom. All the other expectation values are zero essentially because of the orthogonality of spherical harmonics
$$ \int d\Omega\,Y^*_{L'M'} Y_{1M}=0 $$
for $L\neq 1$. Now, the $z$ expectation value is only nonzero for $M=0$ and $x+iy$ and $x-iy$ are only nonzero in the states with $M=\pm 1$, as usual. Here I used that $x,y,z$ may be expressed as combinations of $r Y_{1M}$ for the three values of $M$.
Concerning the remaining factor, the actual result, yes, you may relate it to the expectation value of $r$ (or its power, depending on the conventions for the decomposition into radial and angular functions, and the power of $r$ that you add) but only in the radial part of the wave function. The right formulae are easily obtained by separating the radial and angular variables.
Be careful. The expectation value of $r$ in the whole 3-dimensional wave function is of course zero due to analogous addition-of-angular-momentum arguments: the integral of any $Y_{LM}$, including $L=1$, over the sphere is zero unless $L=0$ because it's the orthogonality again. So only the expectation value of $r$ (or the right power) in the radial part of the wave function is relevant.
