The hyperfine interaction can be written in the form: $${\cal{H}}_{\rm{hf}} = {\mu}_I \cdot \vec{h}_{\rm{eff}}= g_{\rm{N}} \mu_{\rm{N}} \vec{I}\cdot{\vec{h}}_{\rm{eff}}$$ where $\mu_I:=g_{\rm{N}}\mu_{\rm{N}}\vec{I}$ is the nuclear magnetic moment ($g_{\rm{N}}$ being the nuclear g-factor, and $\mu_{\rm{N}}\simeq \mu_B/1836$ the nuclear magneton), and $\vec{h}_{\rm{eff}}$ is the total effective magnetic field produced by the electron. There is an alternative expression that emphasizes the two angular momentum operators involved, nuclear ($\vec{I}$) and electronic ($\vec{j}$), $${\cal{H}}_{\rm{hf}} = a_{\rm hf} \vec{I} \cdot \vec{j}$$ in terms of a coupling hyperfine constant $a_{\rm{hf}}$, which has to be determined.
By direct comparison of these formulas, one deduces $a_{\rm{hf}}\vec{j}=-g_{\mathrm{N}}\mu_{\mathrm{N}}\vec{h}_{\rm{eff}}$. Applying the inner product with $\vec j$, $ a_{\rm{hf}} j^2 =- g_{\mathrm{N}}\mu_{\mathrm{N}}\vec{h}_{\rm{eff}}\cdot\vec{j}$, we see $a_{\rm{hf}}$ can be obtained via the expectation values
$$\displaystyle a_{\rm{hf}} =- g_{\mathrm{N}}\mu_{\mathrm{N}}\frac{\langle \vec{h}_{\rm{eff}}\cdot\vec{j} \rangle}{j(j+1)}$$
on the electronic state denoted by $|j, m_j (l, s)\rangle$. Here's a simple (semi-classical) derivation for the full expression of $\vec{h}_{\rm{eff}}$.
This field has essentially two contributions, one due to the electronic orbital angular momentum $\vec \ell$, and the other given by its spin $\vec s$: $$\vec{h}_{\rm {eff}} = \vec{h}_{\ell} + \vec{h}_s$$ For what concerns $\vec{h}_{\ell}$, the magnetic field perceived by the proton can be written in terms of the electric field $\vec \varepsilon$: $$\vec{h}_{\ell}=\frac{\vec{\varepsilon}\times \vec{v}}{c} = - \frac{e}{r^3} \frac{\vec{r}\times\vec{v}}{c} = - \frac{e}{r^3} \frac{\vec{r}\times m\vec{v}}{m c} $$ Now using the quantization $\displaystyle \vec{r}\times m\vec{v}=:\hbar \vec{\ell}$, and the expression for the Bohr magneton $\displaystyle \mu_B = \frac{e\hbar}{2mc}$ $$\vec{h}_{\ell} = -\frac{2\mu_B}{r^3}\vec{\ell}$$ For $\vec{h}_s$, we can exploit the classical expression for a magnetic dipole: $${\displaystyle \mathbf {B} (\mathbf {m} ,\mathbf {r} )={\frac {\mu _{0}}{4\pi r^{3}}}\left(3(\mathbf {m} \cdot {\hat {\mathbf {r} }}){\hat {\mathbf {r} }}-\mathbf {m} \right)+{\frac {2\mu _{0}}{3}}\mathbf {m} \delta ^{3}(\mathbf {r} )}$$ adapted in this case to: $$\vec{h}_s = \frac{2\mu_B}{r^3}\left(\vec{s}-\frac{3(\vec{s}\cdot \vec{r})\vec{r}}{r^2}\right)-\frac{8\pi}{3}2\mu_B\vec{s}\delta(\vec{r})$$ The term containing the delta function is the so-called Fermi contact interaction when $\vec{r}=0$.
$$\displaystyle \langle \vec{h}_{\rm{eff}} \cdot \vec{j} \rangle = \langle \vec{h}_{\rm{eff}} \cdot \vec{l} \rangle + \langle \vec{h}_{\rm{eff}} \cdot \vec{s} \rangle$$
$$\require{cancel}\displaystyle\langle|\frac{2\mu_B}{r^3}\left(\cancel{\vec{l}\cdot \vec{s}} - l^2 - \frac{3(\vec{s}\cdot\vec{r})(\vec{l} \cdot \vec{r}) }{r^2}+ s^2 - \cancel {\vec{l}\cdot \vec{s}} - \frac{3(\vec{s}\cdot\vec{r})^2} {r^2}\right) - \frac{8\pi}{3}2\mu_B(\vec{s} \cdot \vec{j})\delta(\vec{r})|\rangle$$
Now, $\vec{l}\cdot\vec{r}=0$ (if we reckon the semiclassical expression $\hbar\vec{l}=\vec{r}\times m\vec{v}$), while $\displaystyle \langle(\vec{s}\cdot\vec{r})^2/r^2\rangle = \frac 1 4 \langle (\vec{\sigma}\cdot \hat{r})^2\rangle $. The Pauli vector $\vec{\sigma} \cdot \hat r$ has eigenvalue $\pm|\hat {r}| = \pm 1$ (I'd like to thank @KFGauss for the suggestion), so we end up with:
$$a_{\rm{hf}}= \frac{2\mu_Bg_{\mathrm{N}}\mu_{\mathrm{N}}}{j(j+1)}\langle \frac{l^2}{r^3}+\frac{8\pi}{3}(\vec{s}\cdot\vec{j})\delta(\vec{r})\rangle$$
Finally, here's the question. For $\bf{s}$ states my book says the only term surviving is the Fermi contact interaction. This means $$\langle \frac {l^2}{r^3}\rangle =0$$ But is it?
$$ {\displaystyle\left\langle {\frac {l^2}{r^{3}}}\right\rangle ={\frac {Z^{3}}{n^{3}a_{0}^{3}}}{\frac {\cancel{l(l+1)}}{\cancel{l(l+1)}(l+1/2)}}} = \frac {2Z^{3}}{n^{3}a_{0}^{3}}\neq 0 $$
