Expectation value of $\sum_{a=1}^2\left(\frac{3r_a^ir_a^j}{r_a^5}-\frac{\delta^{ij}}{r_a^3}\right)$ for the helium ground state My question is somewhat related to this one: Expectation value for spherically symmetric states
In the case of the hydrogen atom, we have that in the ground state the following expectation value vanishes,
\begin{equation}
\left\langle\phi(\vec{r})\left|\left(\frac{3r^ir^j}{r^5}-\frac{\delta^{ij}}{r^3}\right)\right|\phi(\vec{r})\right\rangle = 0,
\end{equation}
where $r=|\vec{r}|$.
Does a similar relation also hold for the helium ground state, i.e.
\begin{equation}
\left\langle\phi_\mathrm{He}(\vec{r}_1,\vec{r}_2)\left|\sum_{a=1}^2\left(\frac{3r_a^ir_a^j}{r_a^5}-\frac{\delta^{ij}}{r_a^3}\right)\right|\phi_\mathrm{He}(\vec{r}_1,\vec{r}_2)\right\rangle = 0?
\end{equation}
Here, $\phi_\mathrm{He}(\vec{r}_1,\vec{r}_2)$ is a state with a total angular momentum quantum number $L=0$. However, since it is in general not a product of single-electron wavefunctions the angular momenta of the two individual electrons are not defined (which makes it difficult for me to exploit symmetries for the integration of the two coordinates).
 A: Short story: it's zero, but with interesting physics to explore along the way.


However, since it is in general not a product of single-electron wavefunctions the angular momenta of the two individual electrons are not defined (which makes it difficult for me to exploit symmetries for the integration of the two coordinates).

In principle, yes, it's true that in general there is no guarantee that the wavefunction will be a product of single-electron wavefunctions (or the closest thing for indistinguishable particles, a Slater determinant).
However, the helium ground state is extremely well-described within the Hartree-Fock approximation, which assigns it a ground state of $1\rm s^2$ configuration and $^1\rm S$ symmetry, i.e., with both electrons in a $1\rm s$ orbital. In principle, there are indeed post-Hartree-Fock corrections to this, described either in Configuration Interaction (detailed in e.g. J. Chem. Phys. 30, 617 (1959), Table VIII) or other, fancier methods, and those assign a population of the order of $0.06^2$ to the $2p^2$ term, which is the highest that might be able to contribute to the quadrupole moment.

... but, in any case, the whole thing is moot, and independent of whether the wavefunction is a single Slater determinant or a more correlated state. Ultimately, the ground state has $^1\rm S$ symmetry, which means that it is spherically symmetric, and therefore that any non-symmetric observables must have a zero expectation value.
If you want to make that rigorous, then the tool to use is the Wigner-Eckart theorem, which tells you that the expectation value on an angular-momentum eigenstate $|j,m⟩$ (like the multi-electron $|0,0⟩$ encoded by the $^1\rm S$ symmetry) of a spherical-tensor operator $T_{q}^{(k)}$ (like the quadrupole moment, for which $k=2$) must have the form
$$
\left< j,m \middle| T_{q}^{(k)} \middle| j,m \right>
=
\left< j,m, k,q| j,m \right>
\left( j \middle| T^{(k)} \middle| j \right)
$$
in terms of a reduced matrix element $\left( j \middle| T^{(k)} \middle| j \right)$ that is independent of the orientations $m$ and $q$, and a Clebsch-Gordan coefficient, $\left< j,m, k,q| j,m \right>$, which captures all of the orientation dependence. For your case, this vanishes,
$$
\left< 0,0,2,q| 0,0 \right> = 0
,
$$
as you can't add $k=2$ to $j=0$ and have a final $j=0$, so the matrix element is guaranteed to be zero.
(... but then again, you knew this already.)
