The Qiskit API defines the electronic dipole moment as
$$ \hat{d} = \sum_{p,q} d_{pq}^x a_p^\dagger a_q , $$ where $d_{pq}^x,d_{pq}^y,$ and $d_{pq}^z$ are the Cartesian components of the vector $$ \textbf{d}_{pq} = \int d\tau \phi_p \frac{1}{\textbf{r}} \phi_q . $$
(Link to the API : https://qiskit.org/ecosystem/nature/stubs/qiskit_nature.second_q.properties.ElectronicDipoleMoment.html#qiskit_nature.second_q.properties.ElectronicDipoleMoment)
Can someone please explain how the dipole moment operator can be expressed in the second quantised form this way? I follow that for a discrete one-dimensional system, the expectation value of the dipole moment is easily evaluated by $$ \langle \mu \rangle = \langle \psi | \hat{\mu} | \psi \rangle $$ where $\hat{\mu} = \sum_i \hat{N}_i x_i$ but I cannot extend this to a three-dimensional system with continuous wavefunctions.
Specifically, I cannot see how the term $\frac{1}{\textbf{r}}$ arises in the equation for $\textbf{d}_{pq}$. From dimensional analysis it makes sense to have this form as one needs a factor of $d \tau = 4 \pi r^2 dr$ to convert from Cartesian to spherical polar coordinates, but I'd like to know how it's derived and what happens to the $4 \pi$ factor.